{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import astropy.units as u\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "\n",
    "import warnings\n",
    "from astropy.coordinates import Distance\n",
    "import astropy.constants as const\n",
    "from astropy.time import Time\n",
    "import time\n",
    "from astropy.io import fits\n",
    "import os\n",
    "main_path=os.path.abspath(\"./\")+\"/../\"\n",
    "\n",
    "os.chdir(main_path+\"code/\")\n",
    "from all_script import PlotLC\n",
    "from all_script import LightCurve\n",
    "from all_script import CSP_II_photometry\n",
    "from all_script import match_same_date\n",
    "from all_script import colorLC\n",
    "from all_script import MCMC_Vacca_Leibundgut\n",
    "from all_script import read_ZTF\n",
    "from all_script import spec_color_telescp\n",
    "from all_script import extinction_MW\n",
    "from all_script import Spectra\n",
    "from all_script import curve_fit_mcmc\n",
    "from all_script import MagFlux\n",
    "from all_script import read_autophot_LC\n",
    "from all_script import read_Swift\n",
    "from all_script import read_atlas_0\n",
    "from all_script import superbol_lc2txt\n",
    "import scipy.integrate as integrate\n",
    "from scipy import integrate\n",
    "from all_script import PlotSpec\n",
    "from scipy.optimize import curve_fit\n",
    "from all_script import Airmass"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/ysy/.local/lib/python3.7/site-packages/ipykernel_launcher.py:9: SettingWithCopyWarning: \n",
      "A value is trying to be set on a copy of a slice from a DataFrame\n",
      "\n",
      "See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy\n",
      "  if __name__ == '__main__':\n"
     ]
    }
   ],
   "source": [
    "# https://iopscience.iop.org/article/10.3847/2041-8213/ad527a#apjlad527afn7\n",
    "# 这篇文章下获得的数据。\n",
    "Data_das=pd.read_csv(\"../lc/kaustavkdas/SN2023zaw-main/lightcurve_23bqun_TableC1.csv\")\n",
    "Data_das[\"jd\"]=Time(Data_das[\"jd\"].values,format=\"jd\").mjd\n",
    "Data_das.rename(columns={\"jd\":\"mjd\",\"emag\":\"mag_e\",\"tag\":\"instrument\",\n",
    "                         },inplace=True)\n",
    "replace_filter={'ZTF_g':\"g\", 'ZTF_r':\"r\", 'atlasc':\"c\", 'atlaso':\"o\", 'c':\"c\", 'o':\"o\", 'sdssg':\"g\", 'sdssi':\"i\",\n",
    "       'sdssr':\"r\", 'sdssz':\"z\"}\n",
    "np.unique(Data_das[\"filter\"].values)\n",
    "for f in  replace_filter:\n",
    "    Data_das[\"filter\"][Data_das[\"filter\"]==f]=replace_filter[f]\n",
    "Data_das[\"detection\"]=True\n",
    "Data_das.to_csv(\"../lc/SN2023zaw_lc_Das.csv\",index=None)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAW4AAAD1CAYAAABwdB+7AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAQYklEQVR4nO3df4zkd13H8eebXD1XZXeH3tYWaPtNE1qpGnsyGCupN+UoadCiESy9BEIjOlqJESTGkpDr3v0F1WD5R8umwKHRpkEQWUEJrL2tVsDsebUttNCm7pXSli72MmtMOSV+/GO+e4x3++t2v9+Z+cw+H8lkv/P9fuf7/Xw/O/vaz3zm+/18I6WEJCkfLxp0ASRJ58bglqTMGNySlBmDW5IyY3BLUmYMbknKzK46N75nz55UFEWdu5CkkXPs2LHvpJSm1lpea3AXRcHCwkKdu5CkkRMRJ9ZbbleJJGXG4JakzBjckpQZg1uSMmNwS1JmDG5JyozBLUmZqTW4O50O7Xab2dnZOncjSTtKrcE9MTHBzMwMN9xwQ527GWnFHQVxKJg+Oj3ookgaEnaVbEM/QrWYLNh36T6OPHDEAJcE1HzJ+6grJguKyYLp1nRt+zj5wkkWO4ucP3Y+47vHueaSa2rbl6Q82OLehpMvnOT4s8eZe2Kutn0sdhZZPrXMk50nWT61zMF7D9a2L0l5MLi3YSVU6wzTYqJgfPc4l0xcwvjucQ5fe7i2fUnKg10l21BMFCx2FmsN08ZYg8ZYA4BLJi5h/2X7a9uXpDzY4t6GxliDvRfurTVMW0WL+RPzAMyfmPfLSUlESqm2jTebzTTK43FPH53m0Pwhbtt3W61fUEraWSLiWEqpueZyg1uShstGwW1XiSRlxuCWpMwY3JKUGYNbkjJjcEtSZgxubZsjGEr9lW1wGxbDY2UEQ89ll/oj20ve+zEynyQNo2xb3P0YmU+b4+9C6q9sg7sfI/Npc/xdSP2VbXCvDHfqMKeD5+9C6q9s+7hXhjt1mNPB83ch9Ve2Le6V4U49q0TSTjO0wb3R6X7TrWnSbcmzSoaA/0Sl/hraYV1bR1oAHL35aHUFkqQMbHtY14j4aEQ8FxEP98y7JyIeKB+LEfFAReU9zVPMJGl1m+kqOQJc3zsjpfSWlNJVKaWrgE8Cn6q6YJ5iJkmr2zC4U0r3Ac+vtiwiArgRuLvicnmKmSStYbunA14DfDul9FgVhenlKWaStLrtBvcB1mltLy0t0Wx+v3+93W7Tbre3uUtJ2tm2HNwRsQv4FeBVa60zNTXFls8qKVocmj/E9NFpT/mTpB7baXG/Dng0pfRUVYXpNd0ysCVpNZs5HfBu4EvAFRHxVES8o1x0EzV8KSlJWt+GLe6U0oE15t9ceWkkSRsa2kveJUmrM7glKTMGtyRlxuCWpMwY3JKUGYNbkjJjcEtSZgxujZyVuydNvn9y3bsoSbnK9mbB0lqKyYJisjj93KETNGpscUtSZgxuScqMwS1JmTG4NbK84bRGlcGtkbMS2I8//7g3nNZIMrg1chY7iyyfWobAG05rJHk6oEZOMVGw2FmkmCi84bRGksGtkdMYa9AYawy6GFJt7CqRpMwY3JKUGYNbkjJjcGvktIoW8yfmAZg/Me8gUxo5kVKqbePNZjMtLCzUtn1JGkURcSyl1FxruS1uScqMwS1Jmak1uDudDu12m9nZ2Tp3I0k7Sq0X4ExMTDAzM1PnLiRpx7GrRJIyY3BLUmYMbknKjMEtSZkxuCUpMwa3tI7ijoI4FF42r6HieNzSOorJgmKyYLo1PeiiSKfZ4pakzBjckpQZg1uSMmNwS1JmDG5JyozBLa3j5AsnOf7sceaemBt0UaTTDG5pHYudRZZPLXPw3oODLop0msEtraOYKBjfPc7haw8PuijSaV6AI62jMdagMdZg/2X7B10U6TRb3JKUGYNbkjJjcEtSZgxuScqMwS2to1W0mD8xv+VhXR0WVnWIlFJtG282m2lhYaG27UvDrnWkBcDRm48OtBzKS0QcSyk111pui1uSMrNhcEfERyPiuYh4uGfeVRHx5Yh4ICIWIuJn6i2mJGnFZlrcR4Drz5h3O3AopXQVcLB8Lknqgw2DO6V0H/D8mbOB8XJ6Ani64nJJktaw1Uve3wV8PiL+iG74/1xlJZIkrWurX07eArw7pXQx8G7gI6uttLS0RLPZPP2YmZnZajklSaWttrjfDvxuOf0J4K7VVpqamsLTASWpWlttcT8N7CunXws8Vk1xJEkb2bDFHRF3Ay1gT0Q8BdwG/AbwoYjYBXwXaNdZSEnS920Y3CmlA2sselXFZZEkbYJXTko18p6VqoPBLdXIe1aqDga3VKOq7lnpKIPq5T0npRpVdc/KYrKgmCyYbk1XUzBlzRa3JGXG4JakzBjckpQZg1uSMmNwS1JmDG5JyozBLUmZMbglKTMGt5QBxzxRL4NbqlGraDF/Yn7bl6o75ol6ecm7VKPp1nQll6kXEwWLncVtj3mi0WBwSxmoaswTjQa7SiQpMwa3JGXG4JakzBjckpQZg1uSMmNwS1JmDG5JyozBLUmZqTW4O50O7Xab2dnZOncjSTtKrcE9MTHBzMwMN9xwQ527kUZeVWOeaDRESqm2jTebzbSwsFDb9iVpFEXEsZRSc63l9nFLUmYMbknKjMEtaUuKOwriUNjvPgAO6yppS4rJgmKyqGS8cZ0bW9zSCLI1PNpscUsjyNbwaLPFLUmZMbglKTMGtyRlxuCWtCUnXzjJ8WePM/fE3KCLsuMY3JK2ZLGzyPKpZQ7ee3DQRdlxDG5JW1JMFIzvHufwtYcHXZQdx9MBJW1JY6xBY6zB/sv2D7ooO44tbknKjMEtSZkxuCUpMwa3JGXG4JakzBjckrbE+2AOjqcDStqS6da0ow8OiC1uScqMwS1JmdkwuCPioxHxXEQ83DPvpyLiSxHxUETMRsR4vcWUJK3YTIv7CHD9GfPuAm5NKf0k8NfA71dcLknb4Mh9o23D4E4p3Qc8f8bsK4D7yukvAG+quFyStsGR+0bbVvu4HwbeWE7/KnBxNcWRVAVH7httWw3uXwPeGRHHgBcD/73aSktLSzSbzdOPmZmZrZZT0jlojDXYe+FeR+4bUVs6jzul9CjweoCIuBz4hdXWm5qaYmFhYeulkySdZUst7oi4oPz5IuB9wJ1VFkqStLbNnA54N/Al4IqIeCoi3gEciIhvAI8CTwMfq7eYkqQVG3aVpJQOrLHoQxWXRZLOSXFHwYnOCW7bd9uOuvzesUokZauYLCgmix0V2uAl75KUHYNbkipQ3FEQh6Ivw9zaVSJJFehnt40tbknKjMEtSZkxuCUpMwa3JGXG4JakzBjckpQZg1uSMmNwS1JmDG5pBLWKFvMn5vtyFZ/6zysnpRE03ZrecQMv7SS2uCWpYnWPW2JwS8rWyRdOcvzZ48w9MTfoovw/xWTBvkv31fapx+CWlK3FziLLp5Y5eO/BQRelrwxuSdkqJgrGd49z+NrDgy5KX/nlpKRsNcYaNMYa7L9s/6CL0le2uCUpMwa3JGXG4JakzBjcklSxuk9TNLglqWJ1n6ZYa3B3Oh3a7Tazs7N17kaSBq63lV33aYq1ng44MTHBzMxMnbuQpKHQ28qu+zRFu0okqQL9vBjI4JakCjTGGuy9cG9fLgYyuCUpMwa3JGXG4JakzBjckpQZg1uSKlb3PT8d1lWSKlb3PT9tcUvK1k69m70tbknZ2ql3s7fFLUmZMbglKTMGtyRlxuCWpMwY3JKUGYNbkjJjcEtSZgxuScqMwS1JmTG4JSkzBrckZcbglqTMGNySVIF+jlQYKaX1V4i4GPgz4ELgf4GZlNKHIuIlwD1AASwCN6aUTva+ttlspoWFhRqKLUmjKyKOpZSaay3fTIv7e8B7UkqvBH4WeGdEXAncCsyllF4BzJXPJUk12zC4U0rPpJT+tZz+T+AR4GXALwEfL1f7OPDLNZVRktTjnPq4I6IA9gJfAX40pfQMdMMduKDy0kmSzrLpO+BExI8AnwTelVJajogNX7O0tESz+f1umna7Tbvd3ko5JUmlTQV3RJxHN7T/IqX0qXL2tyPiopTSMxFxEfDcma+bmprCLyclqVobdpVEt2n9EeCRlNIHexZ9Bnh7Of124G+qL15/zMzMDLoIQ8l6OZt1sjrr5Wx11slm+rhfA7wNeG1EPFA+3gC8H7guIh4DriufZ8k33eqsl7NZJ6uzXs5WZ51s2FWSUvonYK0O7f3VFkeStJENL8DZ1sYjloATte2gOnuA7wy6EEPIejmbdbI66+Vs26mTS1NKU2strDW4JUnVc6wSScqMwS1JmRmp4I6IyYj4q4h4NCIeiYirI+IlEfGFiHis/Nko170uIo5FxEPlz9f2bOdAOf/BiPj7iNjTs+zGiPhaRHw1Iv5yEMd5ruqul4i4JCLujYjj5bI3DOpYN6vCOnlLecxfjYjbe+bvjoh7IuLxiPhKedXxUOtDnfxe+bfzYETMRcSlgzjOc1V3vfQsf3NEpIhYc3Cp01JKI/OgO2bKr5fTPwBMArcDt5bzbgU+UE7vBV5aTv8E8K1yehfdi4n2lM9vB6bL6VcAx4FG+fyCQR/zkNTLDHBLOX0lsDjoY+5TnZwPPAlM9Wxzfzn928Cd5fRNwD2DPuYhqJNrgR8qp2/JoU76US/l8xcD9wFfBpoblmnQlVJh5Y4D/075hWvP/K8DF5XTFwFfX+W1AfwHsBs4D1gCLi3n3wm0y/VuX/kF5vLoU718GPiDcvpq4J8Hfdx9qpNXA1/sWfY24E/K6c8DV5fTu+ieXRBVH0tOdXLGa/YC9w/6uIelXoA7gF8Ejm4muEepq+QyusHysfIj+10R8cNsbjCsNwHHU0qnUkr/Q7c18BDwNN0W5EfK9S4HLo+I+yPiyxFxfc3HVIV+1Ms08NaIeAr4HPA7dR5QBSqpE+Bx4MciooiIXXRHyLy4XO9lwDfLbX0P6NBtdQ2rftRJr3cAf1fDcVSt9nqJiL3AxSmlv91soUYpuHcBPw38aUppL/BfbGKM8Ij4ceADwG+Wz8+jG1B7gZcCDwLv7dnHK4AWcAC4KyImqzyIGvSjXg4AR1JKLwfeAPx5RAzze6uSOkndG4fcQveGIv9I94Yi31tZfZVNDPO5t/2ok5XXvBVoAn9YXfFrU2u9lH8nfwy855xKNeiPIhV+pLmQnr5V4Brgs6zzkQZ4OfAN4DU9815N9wYRK89/HvhcOX0ncHPPsjng1YM+9iGol6/SbTGsLHuCIe7/r6pOVtluG7i9nM6tq6T2Oimfv47umP5D+/7oZ70AE+V7Y7F8fJfuJ9p1u0uGuVV0TlJKzwLfjIgryln7ga+xxmBYZUv5s8B7U0r392zqW8CVEbFy1dJ1dN9oAJ+m+wUL5RkVl9MNqaHVp3p5stwuEfFK4AfpfrwcShXWCRFxQfmzQfcLybvKRb3bejPwD6n8ix1G/aiTskvgw8AbU0pnjSY6jOqul5RSJ6W0J6VUpJQKul9OvjGltP6wqoP+j1bxf8ergAW6H+M/DTTo9ivOAY+VP19Srvs+uh97Huh5XFAu+y26ofQgMAucX84P4IPlL+4h4KZBH/OQ1MuVwP3Av5Xrv37Qx9zHOrm7fD98rff9QPef1yfo9m3+C3DZoI95COrki8C3e9b/zKCPeRjq5Yx9HWUTX056ybskZWZkukokaacwuCUpMwa3JGXG4JakzBjckpQZg1uSMmNwS1JmDG5Jysz/Aa8MvoPslz09AAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# 读取数据并画图。\n",
    "SN2023zaw_lc=LightCurve(name=\"SN2023zaw\",\n",
    "                        csv_name=\"../lc/SN2023zaw_lc_Das.csv\",\n",
    "                        redshift=0.0,\n",
    "                        distmod=0,\n",
    "                        ebv_host=[0.0,0.0],\n",
    "                        ebv_MW=[0.0,0.0],\n",
    "                        obs_mag=False)\n",
    "PlotLC(SN2023zaw_lc).plot_band(\"g\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[6.57647502e+044 6.35210689e+044 6.13539348e+044 5.92607364e+044\n",
      " 5.72389511e+044 5.52861427e+044 5.33999577e+044 5.15781233e+044\n",
      " 4.98184441e+044 4.81187994e+044 4.64771411e+044 4.48914909e+044\n",
      " 4.33599379e+044 4.18806367e+044 4.04518044e+044 3.90717192e+044\n",
      " 3.77387181e+044 3.64511947e+044 3.52075975e+044 3.40064278e+044\n",
      " 3.28462381e+044 3.17256304e+044 3.06432541e+044 2.95978051e+044\n",
      " 2.85880234e+044 2.76126922e+044 2.66706362e+044 2.57607201e+044\n",
      " 2.48818474e+044 2.40329590e+044 2.32130319e+044 2.24210781e+044\n",
      " 2.16561432e+044 2.09173054e+044 2.02036744e+044 1.95143902e+044\n",
      " 1.88486221e+044 1.82055679e+044 1.75844526e+044 1.69845278e+044\n",
      " 1.64050704e+044 1.58453823e+044 1.53047890e+044 1.47826389e+044\n",
      " 1.42783029e+044 1.37911733e+044 1.33206629e+044 1.28662048e+044\n",
      " 1.24272514e+044 1.20032737e+044 1.15937607e+044 1.11982189e+044\n",
      " 1.08161718e+044 1.04471589e+044 1.00907355e+044 9.74647215e+043\n",
      " 9.41395394e+043 9.09278018e+043 8.78256384e+043 8.48293108e+043\n",
      " 8.19352082e+043 7.91398430e+043 7.64398468e+043 7.38319656e+043\n",
      " 7.13130570e+043 6.88800854e+043 6.65301190e+043 6.42603258e+043\n",
      " 6.20679707e+043 5.99504117e+043 5.79050969e+043 5.59295617e+043\n",
      " 5.40214254e+043 5.21783885e+043 5.03982301e+043 4.86788050e+043\n",
      " 4.70180411e+043 4.54139371e+043 4.38645600e+043 4.23680427e+043\n",
      " 4.09225816e+043 3.95264351e+043 3.81779206e+043 3.68754130e+043\n",
      " 3.56173428e+043 3.44021939e+043 3.32285020e+043 3.20948526e+043\n",
      " 3.09998797e+043 2.99422638e+043 2.89207303e+043 2.79340482e+043\n",
      " 2.69810285e+043 2.60605228e+043 2.51714218e+043 2.43126541e+043\n",
      " 2.34831847e+043 2.26820141e+043 2.19081769e+043 2.11607406e+043\n",
      " 2.04388043e+043 1.97414982e+043 1.90679820e+043 1.84174439e+043\n",
      " 1.77891001e+043 1.71821933e+043 1.65959923e+043 1.60297906e+043\n",
      " 1.54829058e+043 1.49546790e+043 1.44444735e+043 1.39516747e+043\n",
      " 1.34756885e+043 1.30159415e+043 1.25718795e+043 1.21429675e+043\n",
      " 1.17286886e+043 1.13285436e+043 1.09420502e+043 1.02415087e+043\n",
      " 9.42244398e+042 8.56191885e+042 7.66978893e+042 6.75892601e+042\n",
      " 5.84512376e+042 4.94667419e+042 4.08354057e+042 3.27609385e+042\n",
      " 2.54345514e+042 1.90159186e+042 1.36143534e+042 9.27386810e+041\n",
      " 5.96608314e+041 3.59403998e+041 2.00766589e+041 1.02832771e+041\n",
      " 4.76751917e+040 1.97118639e+040 7.14546643e+039 2.22685120e+039\n",
      " 5.83353832e+038 1.25175492e+038 2.13572735e+037 2.80015384e+036\n",
      " 2.71262018e+035 1.85605403e+034 8.51715127e+032 2.46977585e+031\n",
      " 4.22657109e+029 3.94613234e+027 1.83658872e+025 3.84126840e+022\n",
      " 3.20492330e+019 9.30232665e+015 8.02587095e+011 1.71807092e+007\n",
      " 7.41430744e+001 5.08138818e-005 4.20481653e-012 3.06623682e-020\n",
      " 1.37223658e-029 2.48704173e-040 1.13225267e-052 7.47978573e-067\n",
      " 3.81682024e-083 7.29082443e-102 2.26803088e-123 4.41583773e-148\n",
      " 1.79350802e-176 4.30010988e-209 1.42705537e-246 0.00000000e+000\n",
      " 0.00000000e+000 0.00000000e+000 0.00000000e+000 0.00000000e+000\n",
      " 0.00000000e+000 0.00000000e+000 0.00000000e+000 0.00000000e+000\n",
      " 0.00000000e+000 0.00000000e+000 0.00000000e+000 0.00000000e+000\n",
      " 0.00000000e+000 0.00000000e+000 0.00000000e+000 0.00000000e+000\n",
      " 0.00000000e+000 0.00000000e+000 0.00000000e+000 0.00000000e+000\n",
      " 0.00000000e+000 0.00000000e+000 0.00000000e+000 0.00000000e+000\n",
      " 0.00000000e+000 0.00000000e+000 0.00000000e+000 0.00000000e+000]\n",
      "[0.00000000e+00 2.59531003e+33 5.57717227e+33 9.00315902e+33\n",
      " 1.29394168e+34 1.74619435e+34 2.26580548e+34 2.86280703e+34\n",
      " 3.54872495e+34 4.33680166e+34 5.24225165e+34 6.28255510e+34\n",
      " 7.47779528e+34 8.85104613e+34 1.04288175e+35 1.22415668e+35\n",
      " 1.43242867e+35 1.67171801e+35 1.94664364e+35 2.26251219e+35\n",
      " 2.62542041e+35 3.04237277e+35 3.52141654e+35 4.07179705e+35\n",
      " 4.70413597e+35 5.43063616e+35 6.26531695e+35 7.22428445e+35\n",
      " 8.32604208e+35 9.59184722e+35 1.10461209e+36 1.27169185e+36\n",
      " 1.46364702e+36 1.68418021e+36 1.93754492e+36 2.22862748e+36\n",
      " 2.56304112e+36 2.94723402e+36 3.38861345e+36 3.89568828e+36\n",
      " 4.47823267e+36 5.14747404e+36 5.91630894e+36 6.79955093e+36\n",
      " 7.81421521e+36 8.97984538e+36 1.03188887e+37 1.18571270e+37\n",
      " 1.36241709e+37 1.56540279e+37 1.79857540e+37 2.06642015e+37\n",
      " 2.37408779e+37 2.72749309e+37 3.13342790e+37 3.59969093e+37\n",
      " 4.13523654e+37 4.75034553e+37 5.45682102e+37 6.26821310e+37\n",
      " 7.20007642e+37 8.27026567e+37 9.49927414e+37 1.09106220e+38\n",
      " 1.25313012e+38 1.43922851e+38 1.65291127e+38 1.89825571e+38\n",
      " 2.17993919e+38 2.50332676e+38 2.87457149e+38 3.30072931e+38\n",
      " 3.78989025e+38 4.35132856e+38 4.99567426e+38 5.73510908e+38\n",
      " 6.58359020e+38 7.55710560e+38 8.67396529e+38 9.95513326e+38\n",
      " 1.14246057e+39 1.31098413e+39 1.50422512e+39 1.72577550e+39\n",
      " 1.97974123e+39 2.27081396e+39 2.60435210e+39 2.98647265e+39\n",
      " 3.42415483e+39 3.92535703e+39 4.49914839e+39 5.15585663e+39\n",
      " 5.90723372e+39 6.76664105e+39 7.74925579e+39 8.87230007e+39\n",
      " 1.01552944e+40 1.16203370e+40 1.32924092e+40 1.51997081e+40\n",
      " 1.73740060e+40 1.98510342e+40 2.26708900e+40 2.58784597e+40\n",
      " 2.95238527e+40 3.36628332e+40 3.83572367e+40 4.36753488e+40\n",
      " 4.96922199e+40 5.64898790e+40 6.41574011e+40 7.27907673e+40\n",
      " 8.24924453e+40 9.33705973e+40 1.05537803e+41 1.19109169e+41\n",
      " 1.34199654e+41 1.50920455e+41 1.69374228e+41 1.89648953e+41\n",
      " 2.11810231e+41 2.35891821e+41 2.61884308e+41 2.89721880e+41\n",
      " 3.19267380e+41 3.50296069e+41 3.82478903e+41 4.15366649e+41\n",
      " 4.48376746e+41 4.80785494e+41 5.11728597e+41 5.40212755e+41\n",
      " 5.65139409e+41 5.85339505e+41 5.99618317e+41 6.06815163e+41\n",
      " 6.05892059e+41 5.96067426e+41 5.76994158e+41 5.48946325e+41\n",
      " 5.12944797e+41 4.70746091e+41 4.24655424e+41 3.77192439e+41\n",
      " 3.30702701e+41 2.87033980e+41 2.47367901e+41 2.12230113e+41\n",
      " 1.81632858e+41 1.55268098e+41 1.32678775e+41 1.13373134e+41\n",
      " 9.68838808e+40 8.27915213e+40 7.07304712e+40 6.03883303e+40\n",
      " 5.15022700e+40 4.38537689e+40 3.72623993e+40 3.15793555e+40\n",
      " 2.66813236e+40 2.24651041e+40 1.88431856e+40 1.57402959e+40\n",
      " 1.30908383e+40 1.08370576e+40 8.92775879e+39 7.31741559e+39\n",
      " 5.96553145e+39 4.83615302e+39 3.89746735e+39 3.12144303e+39\n",
      " 2.48349435e+39 1.96216042e+39 1.53879819e+39 1.19729059e+39\n",
      " 9.23772438e+38 7.06375591e+38 5.34994468e+38 4.01071988e+38\n",
      " 2.97405449e+38 2.17971364e+38 1.57768048e+38 1.12674613e+38\n",
      " 7.93250306e+37 5.49959634e+37 3.75071850e+37 2.51335103e+37\n",
      " 1.65272893e+37 1.06506218e+37 6.71654193e+36 4.13848206e+36\n",
      " 2.48736528e+36 1.45568639e+36 8.27936328e+35 4.56708121e+35\n",
      " 2.43804526e+35 1.25657147e+35 6.23712654e+34 2.97349634e+34]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD5CAYAAADcDXXiAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAhSElEQVR4nO3de3iU5Z3/8fd3JifIiZATCQmEMwQ5B1BYK/VXK2oFRZfiqbWlzWrXXt1ft1Xb/e263bbb0253a6WekdrWAx7aYqu1aksRFSWgKBBAQA4xHMIxHARCcv/+yBAmMYEJM5NnZvJ5XRfXleeemWe+N5P58PCde57HnHOIiEhi8XldgIiIRJ7CXUQkASncRUQSkMJdRCQBKdxFRBKQwl1EJAEleV3AKXl5ea6srMzrMkRE4sqKFSv2OOfy247HTLiXlZVRVVXldRkiInHFzLa2N662jIhIAopKuJtZupmtMLPPBLanmdmrZnafmU2LxnOKiMhpIYW7mc03s91mtrrN+HQzW29mG83szqCb7gAWBm074DCQBtSEW7SIiJxZqEfuC4DpwQNm5gfmAZcB5cB1ZlZuZp8C1gK7gu7+qnPuMppD/zvhFi0iImcW0geqzrklZlbWZngSsNE5txnAzJ4AZgIZQDrNgf+RmT3vnGsKPGY/kNrec9TV1VFRUdGyXVlZSWVlZSemIiIip4SzWqYvsD1ouwaY7Jy7DcDMbgb2OOeazGwWcCnQC7invZ3l5+drtYyISISEE+7WzljL+YOdcwuCfn4WeDaM5xIRkU4IJ9xrgNKg7RKgNrxyOu+hVzfzp9U7KS/Oorwoi/LiLIYWZpKW7O/qUkREYkY44b4cGGJmA4APgTnA9RGpqhOy0pIxg2dXfsijx5vX8vt9xqD89JawLy/Kprw4i97pKV1dnoiIJ0IKdzN7HJgG5JlZDXCXc+5hM7sNeBHwA/Odc2uiVmkHZk8sZfbEUpqaHNv3H2VtbT1rd9SztraeZZv38bt3Tv9nok9WWssR/sji5uAvzemJz9deh0lEJH5ZrFxmr6KiwkXjA9W9h49TveMQa3ccbAn+TXVHaGxqnndGahIjijJbHeUPKcxQW0dE4oKZrXDOVbQdj5lzy0RLbkYqfzcklb8bktcydqyhkQ27DrU6yn96RQ1H3mgEmts6g/MzWvXxy4uyyFFbR0TiRMKHe3vSkv2MLunF6JJeLWNNTY5t+462hP3aHfW8sWkvv337w5b7FGWntYT9yMBRfklOD7V1RCTmdMtwb4/PZ5TlpVOWl87lo4paxvccPk51UOCvra3nr+t3E+jqkJmaxIigo/vy4iyGFGaQmqS2joh4R+F+FnkZqVw4JJ8Lh5w+XfKxhkbW7zzU6ih/YdV2jp5obusk+YzBBRmtWjrlxVn06qm2joh0DYX7OUhL9jOmtBdjSnu1jDU1ObbuO7Vap/nD29c27eHZoLZOcXZgtU5xdsuKnZKcHpiprSMikaVwjxCfzxiQl86AvHSuGN26rRPc0lm7o56/rGvT1mnzwa3aOiISLoV7lOVlpPKJofl8Yujpts5HJxpZ37Jap/ko/8nl2/mooU1bp03oq60jIqFSuHugR4qfsaW9GBvU1mlscmzde6TVEf7S9/fw7MrTbZ2+vXq0fHg7MhD4auuISHsU7jHC7zMG5mcwMD+Dz4wubhmvO9S8WmdNS2vnIK+s28Wp755lpiV97IPbIQWZpCTpCooi3ZnCPcblZ6aSn9m6rXP0xMmPrdZ54q3TbZ1kvzG4ILN16Bdlkd0z2atpiEgXU7jHoZ4pSYzrl8O4fjktY41Nji17j7T68HbJ+3U8s/L0VQ379urxsXPr9O2lto5IIlK4J4jmM2FmMCg/gyvHnG7r7D50jOodh1hTe/rcOi9Xn27rZKUltTpzZnlRFoMLMtTWEYlzCvcEV5CZRkFmGhe1aeus29n63DqPvbWVYw3NV0NM9htDCjJb9fFHFGWR3UNtHZF4oXDvhnqmJDG+Xw7j27R1PtjTerXO4vW7eXrF6bZOSU6PoHPrNB/pF2enqa0jEoMU7gIEzoRZkMHgggxmtGnrnAr7NbX1VNfW81JQWye7R/LHVusMLsgg2a+2joiXFO5yRgWZaRQMS2PasIKWsSPHA22doKP8Xy/byvGTzW2dFL+PIYWtz60zojiLrDS1dUS6isJdOi09NYkJ/XOY0P90W+dkYxNb9h4JWo/ffJqFp4LaOqW9A22douyW1TpFauuIRIXCXSIiye9jcEEmgwsymTm2LwDOOeoOHWfNqSP8QPD/ee3ptk6vnskt6/DLA4E/KF9tHZFwKdwlasyMgqw0CrLS+GRQW+fw8ZOs39n6HPm/atPWGdonIyj0sxlRlEmm2joiIVO4S5fLSE1iQv/eTOjfu2XsZGPTx1brvFy9m4VVp9s6/Xr3bNXHH9k3iz5ZauuItEfhLjEhye9jSGEmQwpbt3V2HzoetFqn+YtYf1qzs+VxOT2T25w9M5uB+elq60i3p3CXmGVmFGalUZiVxieHt27rrNvR+hz5v3xjKydOtXWSfAwrzGzVxx/eR20d6V4U7hJ3MlKTqCjrTUVZ67bO5j2tz63z57U7ebJqe8t9+uf2bOnjj+zbfJRfmJWqto4kJIW7JIQkv4+hhZkMLczkqnGn2zq76o+3XBDlVOi/sPp0W6d3esrHvoQ1MC+dJLV1JM4p3CVhmRl9stPok53GxcMLW8YPHWs4fW6dQOgveH1Lq7bO8D6tT5k8vCiLjFS9XSR+mDu14NhjFRUVrqqqyusypJtqaGxic92RVkf5a2rrOXC0oeU+Zbk9g06Z3HxunYJMtXXEW2a2wjlX0XZchyIiQLLfx7A+mQzrk8nV45rHnHPsrD/W6gh/TW09z793uq2Tm57ysWvdDlBbR2KAwl2kA2ZGUXYPirJ78H9GnG7r1B9rYN2OQ6ytPdjcx99RzyOvbeFEY3NbJ/VUWyco9If3ySJdbR3pQlFpy5hZOrAEuMs59wczGwF8DcgDXnHO3dv2MWrLSDxraGxiU93hjx3lH/youa1jBmW56a2/hFWcRb7aOhKmsNoyZjYf+Ayw2zl3XtD4dOBngB94yDn3w8BNdwALT93POVcN3GJmPuDBc56FSIxK9vsY3qf5CH3W+OYx5xw7Dh5rtVLnvQ8P8sf3drQ8Li8jhRFtrnWrto5EQqj/T1wA3AM8emrAzPzAPOASoAZYbmaLgGJgLZAWvAMzmwHcGdiPSMIzM4p79aC4Vw8+Vd66rVNd2/pLWPOXfkBDY/P/olOSfEwdlMtV4/py+agifdtWzknIbRkzKwP+cOrI3cwuAP7dOXdpYPtbgbtmAOlAOfARcLVzriloP390zl3Rdv9qy0h3duJkc1unekc9qz+s58U1O/nwwEcUZ6cx98KBzJlYqp69tKujtkw44X4tMN0596XA9k3AZOfcbYHtm4E9gZ77NGAWkAq865yb13b//fv3d/n5p6/zWVlZSWVlZSemKJI4mpocizfs5r7Fm3lryz5yeiZzx/ThzK4oxedTj15Oi8ZSyPZ+w1r+pXDOLQj6eTGw+Ew7y8/PR0fuIs18PuPi4YVcPLyQFVv388MXqrnz2fd4Yvl2vjvzPEaVZHtdosS4cJp5NUBp0HYJUBteOSLS1oT+OSz8hwv46ewx1Oz/iJnzlnL3K+/T2BQbX0CU2BROuC8HhpjZADNLAeYAiyJTlogEMzNmjS/hlX++iBljivnpSxu4+ZG3OHD0hNelSYwKKdzN7HHgDWCYmdWY2Vzn3EngNuBFoBpY6JxbE71SRSS7RzL/89mx/OiaUSzbvJer5r3G1r1HvC5LYpDOLSMSp6q27ONLj1aR7Pfxq7mTGN4ny+uSxAMdfaCqBbQicaqirDdP/cMF+M24/sE3eX/XIa9LkhiicBeJY0MKM3m88nySfMb1D73J9n1HvS5JYoTCXSTODchL5zdfmszxhka+uGA59ccazv4gSXgKd5EEMKQwk3tvnMAHe47w1cfe5mRj09kfJAlN4S6SIKYOzuM/Zp7H3zbU8b0/VntdjnhMJ6sQSSDXT+7HprrDPLz0A8qLspg9sfTsD5KEpCN3kQTz7ctHMHVwLv/6+9WsqT3odTniEYW7SILx+4yfzRlHTs8Ubv31ypYLhkj3onAXSUB5GanMu2EctQc+4htPrSJWvqwoXUfhLpKgJvTvzbcuH8FLa3dx/5LNXpcjXUzhLpLAvji1jCtGFfHjP61j2ea9XpcjXUjhLpLAzIwfXjOKstx0vvKblWzZo5OMdRcKd5EEl5mWzEOfr8A5x+cfeYvdh455XZJ0AYW7SDcwMD+Dh2+eyO7641x77xt8oCP4hKdwF+kmxvfL4fHK8zl8/CQz71nKwqrtWkWTwBTuIt3I2NJe/PYrUxjeJ4vbn36Xy372Kk+8tU1XdEpAuliHSDfU1OR49u0PeXDJZtbvOkSSz5gyOI/LzuvDp8sLyc1I9bpECVFHF+tQuIt0Y8453vvwIM+/t5MXVu9g696j+AzOH5jLFaOLmDWuhB4pfq/LlDNQuIvIGTnnqN5xiBdW7+D593awqe4IvdNT+Mq0QXx+ShnJfnVxY5HCXUQ6pWrLPu7+y0aWbKhjSEEG/zHzPC4YlOt1WdKGrqEqIp1SUdabR784iQc/V8HRE41c9+Aybn96FUdPnPS6NAmBwl1EzuiS8kJe/vpF3DptEE+tqGHGPa+xbme912XJWSjcReSseqT4uWP6cH49dzIHjjYw857XWLh8u9dlyRko3EUkZFMH5/HC1y5kYllvbn/mXb7z3BpdrzVGKdxFpFPyM1NZ8IWJfHHqAB55bQtfWLBcFwSJQQp3Eem0JL+Pf7uynB9dM4plm/dy9S9eY3PdYa/LkiAKdxE5Z5+d2I/ffOl8Dhxt4Kp5r/H6xj1elyQBCncRCcukAb35/T9OpU92Gp+b/5Y+aI0RUQl3M0s3sxVm9pnA9kAze9jMno7G84mIt0p79+TpW6dwwaBcbn/mXX7y4jqammLjC5LdVUjhbmbzzWy3ma1uMz7dzNab2UYzuzPopjuAhac2nHObnXNzI1OyiMSirLRk5t88kesmlTLvr5v42pPvcKyh0euyuq1Qj9wXANODB8zMD8wDLgPKgevMrNzMPgWsBXZFsE4RiQPJfh//efUo7rxsOM+tquWGh95k7+HjXpfVLYUU7s65JcC+NsOTgI2Bo/ITwBPATOCTwPnA9cCXzUx9fZFuxMy45aJB/OKG8az+8CCz7n2dTVpJ0+XCCd6+QPAnJzVAX+fcvzjn/gl4DHjQOddkZrlmdh8wzsy+1d7O6urqqKioaPnzwAMPhFGaiHjt8lFFzVd+OnaSWb94nWWb93pdUreSFMZjrZ2xlk9QnHMLgn7eC9xypp3l5+ejs0KKJJbx/XL43T9O5eZH3uKmh9/kR9eMZtb4Eq/L6hbCOXKvAUqDtkuA2vDKEZFEU9q7J8/eOpWJZb35+sJV/PSlDbp2axcIJ9yXA0PMbICZpQBzgEWRKUtEEkl2z2QWfGESfz+hhLtfeZ//++Q7HD+plTTRFOpSyMeBN4BhZlZjZnOdcyeB24AXgWpgoXNuTfRKFZF4lpLk48fXjuablw7jd+/UctNDb7H/iC7MHS26EpOIdLlFq2r5xlOr6NurB4/cPJGyvHSvS4pbuhKTiMSMGWOKeexLkzlw9ARX/+I1lm9pu9JawqVwFxFPVJT15rdfmUpOzxRuePBNFq3SeoxIUriLiGfK8tJ59itTGNuvF1974m1++foWr0tKGAp3EfFUr54pPPrFSXxqRCF3LVqjpZIRonAXEc+lJfu594bxzK5oXir5/363mkadVTIs4XxDVUQkYpL8Pn50zWhyM1K5d/Em9h89wf98diypSX6vS4tLCncRiRlmxh3Th5ObnsL3/ljNsYaV/OKG8aQlK+A7S20ZEYk5X7pwIP959Sj+sm43X360SueFPwcKdxGJSddP7sePrx3N0o17+OKC5Rw9cdLrkuKKwl1EYtbsilJ+OnsMyzbv5eb5yzl8XAEfKoW7iMS0q8eV8LM541ixbT83z39LR/AhUriLSMy7ckwxP79uHCu37afy0RXqwYdA4S4iceHyUUX8+NoxLN24h9see5uGxiavS4ppCncRiRvXTijhuzNH8nL1Lr6+cJW+6HQGWucuInHlpgvKOHKikR++sI6eyX5+MGsUPl97V/3s3hTuIhJ3brloEEeOn+Tnf9lITnoKd1423OuSYo7CXUTi0tcvGcq+Iye472+b6JvTg5vO7+91STFF4S4iccnM+M6Mkew8eIy7fr+aPllpXFJe6HVZMUMfqIpI3Ery+/j59eM4r282X318Jau2H/C6pJihcBeRuNYzJYmHPz+R/MxU5v5yOdv2HvW6pJigcBeRuJefmcqCL0ziZJPj5kfe4uBHDV6X5DmFu4gkhEH5Gdx/4wS27TvK1554u9uvgVe4i0jCmDwwl7tmjGTx+jp++tJ6r8vxlFbLiEhCuXFyP9bWHmTeXzdRXpTNFaOLvC7JEzpyF5GEYmb8+4yRjO/Xi288tYrqHfVel+QJhbuIJJzUJD/33TiBzLQkKn9Vxf4jJ7wuqcsp3EUkIRVkpXH/TRPYdfA4//TkOzR1sw9YFe4ikrDG9cvhX68s528b6rh/yWavy+lSEQ93MxthZveZ2dNmdmtgrNzMFprZvWZ2baSfU0SkIzdO7scVo4r4rz+vZ8XWfV6X02VCCnczm29mu81sdZvx6Wa23sw2mtmdAM65aufcLcBsoCJw18uAnzvnbgU+F8H6RUTOyMz4wTWj6NurB1997O1u038P9ch9ATA9eMDM/MA8moO7HLjOzMoDt80AlgKvBO7+K2COmf0EyA2/bBGR0GWlJTPv+vHsOXyC2595F+cSv/8eUrg755YAbf8/MwnY6Jzb7Jw7ATwBzAzcf5FzbgpwQ2B7t3PuH4E7gT2RKl5EJFSjSrL55qXDeGntLhZWbfe6nKgL50tMfYHgv6EaYLKZTQNmAanA8wBmVgZ8G0gHftLezurq6qioqGjZrqyspLKyMozyRERam/t3A/jLut1857m1nD8wl/656V6XFDXhhHt717VyzrnFwOI2g1uAMyZ1fn4+VVVVYZQjInJmPp/x37PHcOn/LuHrC1fxZOX5JPkTc9FgOLOqAUqDtkuA2vDKERGJruJePfjeVeexYuv+hF4eGU64LweGmNkAM0sB5gCLIlOWiEj0zBzblytGF/Gzl9/n/V2HvC4nKkJdCvk48AYwzMxqzGyuc+4kcBvwIlANLHTOrYleqSIikfOdGSNJT/Vz+zPvJuTpgUNdLXOdc67IOZfsnCtxzj0cGH/eOTfUOTfIOff96JYqIhI5eRmp3HXlSN7edoAFr2/xupyIS8xPEkREQjBzbDEXDy/gv15cn3CX51O4i0i3ZWZ8/+rz8PuMb//2vYT6cpPCXUS6taLsHnzz0mEs3biHP763w+tyIkbhLiLd3o3n92dkcRbf/cNaDh8/6XU5EaFwF5Fuz+8zvnvVeeyqP87dr7zvdTkRoXAXEQHG98thzsRS5i/9gA0JsPZd4S4iEnD79OFkpCXxb79fHfcfrircRUQCeqen8M+XDGXZ5n28tHaX1+WEReEuIhLkukn9GJSfzg9fWEdDY5PX5ZwzhbuISJAkv49vXz6CzXuO8Nib27wu55wp3EVE2rh4eAFTBuXyvy9v4OBHDV6Xc04U7iIibZgZ/3LFCA581MAv/rrR63LOicJdRKQdI4uzmTWuhEde38LOg8e8LqfTFO4iIh34p08NwTnHPX+Nvy82KdxFRDpQ2rsnn51YypPLt7N9X3ydNVLhLiJyBrd9cghmxs//El9H7wp3EZEz6JOdxo2T+/PMyg/5YM8Rr8sJmcJdROQsbp02iBS/j5+9vMHrUkKmcBcROYv8zFQ+d0F/Fq2qZUucHL0r3EVEQjD3wgEk+X3cv2Sz16WEROEuIhKCgsw0/n5CCc+sqGF3feyve1e4i4iEqPITAznZ1MTDSz/wupSzUriLiISof246nxldzK+XbeXg0dg+54zCXUSkE265aBBHTjTyq2VbvC7ljBTuIiKdUF6cxSeH5fPIa1s41tDodTkdUriLiHTSly8cyN4jJ3huVa3XpXRI4S4i0kkXDMplaGEGC17fErPXWlW4i4h0kplx85QBrKmtp2rrfq/LaVfEw93MRpjZfWb2tJndGhi7MDD2kJm9HunnFBHpaleNKya7RzILXtvidSntCinczWy+me02s9Vtxqeb2Xoz22hmdwI456qdc7cAs4GKwNirgbE/AL+M7BRERLpez5Qk5kws5U9rdlJ74COvy/mYUI/cFwDTgwfMzA/MAy4DyoHrzKw8cNsMYCnwSpv9XA88Hka9IiIx48bz++Oc49fLtnpdyseEFO7OuSXAvjbDk4CNzrnNzrkTwBPAzMD9FznnpgA3nLqzmfUDDjrn6iNSuYiIx0p79+SS8kIef2tbzC2LDKfn3hfYHrRdA/Q1s2lmdreZ3Q88H3T7XOCRjnZWV1dHRUVFy58HHnggjNJERLrGTeeXsf9oA39eu8vrUlpJCuOx1s6Yc84tBha3c8NdZ9pZfn4+VVVVYZQjItL1pgzKpSSnBwuXb2fGmGKvy2kRzpF7DVAatF0CxO6KfhGRKPD5jNkVpSzduCemrrMaTrgvB4aY2QAzSwHmAIsiU5aISPy4dkIJPoOFVdvPfucuEupSyMeBN4BhZlZjZnOdcyeB24AXgWpgoXNuTfRKFRGJTcW9enDR0HyeqqrhZGOT1+UAoa+Wuc45V+ScS3bOlTjnHg6MP++cG+qcG+Sc+350SxURiV2fnVjKzvpjLHm/zutSAJ1+QEQkIi4eXkheRgpPLo+N1ozCXUQkAlKSfFwzvoRXqnez5/Bxr8tRuIuIRMqs8SWcbHI8/94Or0tRuIuIRMqwPpkM75PJ797+0OtSFO4iIpE0c2xfVm47wLa93q55V7iLiETQlWOKAHjuXW+/06lwFxGJoJKcnkwsy+F3b3/o6VWaFO4iIhE2Y2xf3t99mOodhzyrQeEuIhJhV4wqIsln/H6Vdx+sKtxFRCKsd3oKnxiaz3Pv1NLU5E1rRuEuIhIFM8YUU3vwGCu3eXMBbYW7iEgUXDyigGS/8eKanZ48v8JdRCQKstKSmTo4jxfX7PJk1YzCXUQkSi4d2Ydt+46ybmfXr5pRuIuIRMkl5YWYwZ9Wd31rRuEuIhIleRmpTOzf25O+u8JdRCSKLj2vD+t2HmLr3iNd+rwKdxGRKPp0eSFAlx+9K9xFRKKotHdPzuub1eV9d4W7iEiUXVreh5XbDrC7/liXPafCXUQkyi4Z2dyaWby+6y6erXAXEYmyYYWZ9MlKY/GG3V32nAp3EZEoMzOmDcvn1Q17aGhs6pLnVLiLiHSBacPyOXT8JCu3ds2JxBTuIiJdYOrgPJJ8xuINXdN3V7iLiHSBzLRkJvTP6bIPVRXuIiJdZNqwAqp31LOrC5ZEKtxFRLrItGH5APytC47eFe4iIl1keJ+uWxIZ8XA3sxFmdp+ZPW1mtwbGppnZq4HxaZF+ThGReGBmXDS0a5ZEhhTuZjbfzHab2eo249PNbL2ZbTSzOwGcc9XOuVuA2UBF4K4OOAykATWRK19EJL6cWhL59rYDUX2eUI/cFwDTgwfMzA/MAy4DyoHrzKw8cNsMYCnwSuDurzrnLgPuAL4TftkiIvFpyqA8zOD1TXui+jwhhbtzbgmwr83wJGCjc26zc+4E8AQwM3D/Rc65KcANge1T///YD6RGonARkXiU3TOZkcVZvLFpb1SfJymMx/YFtgdt1wCTAz31WTSH+PMAZjYLuBToBdzT3s7q6uqoqKho2a6srKSysjKM8kREYtOUQXkseG0LxxoaSUv2R+U5wgl3a2fMOecWA4vbDD4LPHumneXn51NVVRVGOSIi8eGCgbk8sGQzK7buZ+rgvKg8RzirZWqA0qDtEqA2vHJERBLfxAG98fssqn33cMJ9OTDEzAaYWQowB1gUmbJERBJXRmoSo0uyo9p3D3Up5OPAG8AwM6sxs7nOuZPAbcCLQDWw0Dm3JmqViogkkCmDcllVc5DDx09GZf+hrpa5zjlX5JxLds6VOOceDow/75wb6pwb5Jz7flQqFBFJQBcMzKOxybF8S9uFiJGh0w+IiHhgQv8cUvw+lkWpNaNwFxHxQI8UP2P79eJ1hbuISGK5YGAua2oPcvBoQ8T3rXAXEfHIlEG5NDl484PIH70r3EVEPDK2Xy++d9V5jCntFfF9h/MNVRERCUNqkp8bz+8flX3ryF1EJAEp3EVEEpDCXUQkASncRUQSkMJdRCQBJUS4P/DAA16XEDGaS+xJlHmA5hKrojGXmA335557LuSfw/mLCd7fudynvdvajsXDXDo7j7bbp34OHovHuUT6NTlTnaHcJ1F+vzq6LR7nEmvvlY4kRLhH6nnO5T6R/oUNR6z9woYjlsI9XLEU7uHQe6XjcS/fKx0x51xUnyBUZlYHbA0aygYOhvhzHnCulzQJ3t+53Ke929qOxcNcOjuPttunfg4ei8e5RPo1OVOdodwnUX6/OrotHucSa++V/s65/LaDMRPuIiISOTHblhERkXOncBcRSUAKdxGRBJTw4W5mI8zsPjN72sxu9bqecJjZVWb2oJn93sw+7XU958rMBprZw2b2tNe1nAszSzezXwZeixu8ricc8f5aBEug90dkMss5F7N/gPnAbmB1m/HpwHpgI3BniPvyAQ8nyFxyvJpLhOfxtNe/Y+cyL+Am4MrAz096XXskXqNYei0iMBfP3h8RnkdYmeX5pM/yF/IJYHzwXwjgBzYBA4EUYBVQDowC/tDmT0HgMTOA14Hr430ugcf9NzA+AeYRM4HSyXl9CxgbuM9jXtcezlxi8bWIwFw8e39Eah6RyKyYvliHc26JmZW1GZ4EbHTObQYwsyeAmc65HwCf6WA/i4BFZvZH4LEoltyhSMzFzAz4IfCCc25llEtuV6Rek1jTmXkBNUAJ8A4x2Nrs5FzWdnF5ndKZuZhZNR6/PzrS2dckEpkVc7+YIegLbA/argmMtcvMppnZ3WZ2P/B8tIvrpE7NBfgq8CngWjO7JZqFdVJnX5NcM7sPGGdm34p2cWHoaF7PAteY2b1AdL9mGDntziWOXotgHb0usfr+6EhHr0lEMiumj9w7YO2MdfhNLOfcYmBxtIoJU2fncjdwd/TKOWedncdeIB7efO3Oyzl3BPhCVxcTpo7mEi+vRbCO5hKr74+OdDSPxUQgs+LxyL0GKA3aLgFqPaolXIkyl0SZR1uJNC/NJfZEdR7xGO7LgSFmNsDMUoA5wCKPazpXiTKXRJlHW4k0L80l9kR3Hl5/inyWT5gfB3YADTT/Kzc3MH45sIHmT5r/xes6u9NcEmUeiTwvzSX2/ngxD504TEQkAcVjW0ZERM5C4S4ikoAU7iIiCUjhLiKSgBTuIiIJSOEuIpKAFO4iIglI4S4ikoAU7iIiCej/A+4TRASr9ZQOAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXMAAAD5CAYAAADV5tWYAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAciElEQVR4nO3dd3hUZd7G8e8zk4QQegkshF4lUgRCjywoKrhBBBFBwA7SbWt/93XXfe2r6xJqVhAVBRFhIdiwYAlNCL0TUDQiJBTpEEjO+0eim2UTmDCTnDNn7s915brISWbm9+wkt9lz7nnGWJaFiIgEN4/dA4iIiP8U5iIiLqAwFxFxAYW5iIgLKMxFRFxAYS4i4gJhdj1w1apVrXr16tn18CIiQSc1NfWAZVnRBX3NtjCvV68eq1evtuvhRUSCjjFmT2Ff02kWEREXUJiLiLiAwlxExAUU5iIiLqAwFxFxAYW5iIgLBF2Yn8vOsXsEERHHCaowP5WVzU2TlzE95Tu0D7uIyL8FPMyNMR5jzDPGmERjzO2BvO9zOTlUKx/J04u2MHbWWk6cORfIuxcRCVo+hbkxZroxJsMYs+m84z2NMduNMWnGmMfyDvcBYoCzQHoghy0XGc7UIW15pGdTPtz4MzdOXEpaxvFAPoSISFDy9S/zGUDP/AeMMV5gItALiAUGGWNigabAcsuyHgRGBm7UXB6PYVS3Rrx1dwcOnciiz4QUPtr4c6AfRkQkqPgU5pZlfQ0cOu9weyDNsqzdlmVlAbPJ/as8HTic9z3ZgRr0fF0aVWXRuHgaVy/HyLfX8OyHW3VxVERClj/nzGOAH/N9np53bB5wnTEmEfi6sBtnZmYSFxf320dSUlKRB6hRoTTv3tuRoR3rkvT1bga/tpKMY6eLfD8iIsHOn10TTQHHLMuyTgJ3X+zG0dHRAdk1sVSYl7/e2Jw2dSvy+LyNJIxPYdLgNsTVq+z3fYuIBAt//jJPB2rn+7wWsNe/cS5d39a1mD+qC1ERXgYmrVB9UURCij9hvgpobIypb4yJAAYCC3298ZEjRxg+fDjJycl+jPCfmtUoz4Ix8XRrWo2nF21h3Ox1qi+KSEgwvvz1aoyZBXQDqgL7gacsy5pmjLkeeBXwAtMty3rG1weOi4uziuvNKXJyLCZ/tYuXF2+nYXRZpgxtS8PossXyWCIiJcUYk2pZVlyBX7PrVERxhvmvUnYeYNzstWSdy+Gl/i3p1aJGsT6eiEhxulCYB9XL+YsqvnFVFo2Np2G1sox8ew3Pqb4oIi5lW5gXxznzgtSsWJo593ZkSMc6TP16N0OmrSTz2JlifUwRkZLm6tMs55u3Jp0n5m+kQulwJg1uQ9u6qi+KSPAI2dMs5+vXphbzRnYhMtzLLVNXMGOp6osi4g4hFeYAsTXLs3BMPN2aRvPn5C3cN3sdJ7NUXxSR4BZyYQ5QoXQ4SUPjePi6pizasJcbJy5ld6Z2XxSR4OX6C6CF8XgMo7s34s27OnDgeBY3TFjKx5v22TKLiIi/QuoCaGF++uUUo2amsj79CPf+vgEPX9uUMG9I/p8WEXEwXQC9iJiKpZkzohODO9Rh6le7GTrtW9UXRSSoKMzzlArz8kzfFrx8cyvW/HCYhMRvSN1z+OI3FBFxAIX5eW5qW4t5ozpTKszLwKTlvLHse9UXRcTxQvYC6IVcXrMCyWPi6do4mqcWbub+d1VfFBFn0wXQC8jJsZi4JI1XPttBk2rlmDykDQ20+6KI2EQXQC+Rx2MYe3Vj3rizPRnHTtNnwlI+2az6oog4j8LcB12bRJM8Np760WW4961Unv9om3ZfFBFHUZj7qFalKN4b0YlbO9Rhyle7uG36txw4rvqiiDiDwrwISoV5ebZvC17q35LUPYdJGJ/Cmh9UXxQR+6nNcglujqvNvFGdCQ8z3DJ1OW8uV31RROylNosfjpw8ywNz1vHFtgz6to7hmb7NiYoIs3ssEXEptVmKSYWocF67LY6HrmnCv9b9RL9Jy/juwAm7xxKREKQw99Ov9cUZd7Zn39HT3JCYwmLVF0WkhCnMA+T3TaJZlFdfHP5WKi98rPqiiJQchXkA1aoUxZx7OzGofR0mf6n6ooiUHIV5gEWGe3muXwte7N+S1XsO0zsxhbWqL4pIMVM1sZgMiKvNvJGd8XoMA6Yu5y3VF0WkGKmaWMx+OZnFA++uY8n2TPq1juGZvi0oHeG1eywRCUKqJtqoYlQE025vxwM9mjB/3U/0nbSU71VfFJEAU5iXAI/HcF+Pxrx+Rzv2HT1N7wkpfLplv91jiYiLKMxLULem1UgeE0+9KmUY9uZqXvpkG9k5Oo8uIv5TmJew2pVzd18c2K42E5fs4vbp33JQ9UUR8ZPC3AaR4V6ev6klL97Ukm+/P0TvxBTW/fiL3WOJSBBTmNtoQLvc+qLHY7h5yjJmrtij+qKIXBKFuc2ax1Rg0dh4ujSqyv/8axMPvbeeU1nZdo8lIkFGYe4AFaMimH57O+7v0Zj5a3Pri3sOqr4oIr7TK0AdwuMx3N+jCa/f0Y6fj5wmITGFz1RfFBEf6RWgDvTjoZOMfDuVTT8dZUz3RjxwTRO8HmP3WCJiM70CNMjUrhzF3BGduSWuNhOWpHHH699y6ESW3WOJiIMpzB0qMtzLC/1b8ny/Fqz87hAJ479RfVFECqUwd7iB7evw/ojOGGMYMGU5b69UfVFE/pvCPAi0qJVbX+zUsApPzt/EH9/bwOmzqi+KyL8pzINEpTIRTL+jHfdd3Zh5a9PpN2kZPxw8afdYIuIQCvMg4vUYHrimCdNvb8dPv5wiIfEbPt+q+qKIKMyDUvfLqrFobDy1K0dx9xureXnxdu2+KBLiFOZBqnblKN4f2ZkBcbVI/EL1RZFQpzAPYpHhXl7s3+q3+mLvxBTWq74oEpIU5i4wsH0d5o7oBMDNU5bzzsofVF8UCTEKc5doWasii8bG07FhFZ6Yv5GH56q+KBJKFOYuUqlMBK/f0Y5xVzdmbqrqiyKhRLsmuozXY3jwmtzdF9MPnyQh8Ru+2Kb6oojbaddEF/vh4ElGzExly89HGXdVI+7rod0XRYKZdk0MUXWqRDFvVGf6t63F+C/SuHPGKg6rvijiSgpzl4sM9/JS/5Y827cFK3YdJCExhQ3pv9g9logEmMI8BBhjuLVDHd7Lqy/2n7yc2d/+YPNUIhJICvMQ0qp2RZLHxtOhQWUem7eRR+auV31RxCUU5iGmcpkIZtzZnrFXNWLO6nRumryMHw+pvigS7BTmIcjrMTx0bVOm3R7Hj4dOkpCYwpJtGXaPJSJ+UJiHsKubVWfR2CupWbE0d72xilc+3aHdF0WClMI8xNWpEsX8UZ3p17oW4z/fyV2qL4oEJYW5EBnu5W83t+SZvs1Znldf3Jh+xO6xRKQIFOYC5NYXB3eoy3sjOmFZFjdNWab6okgQUZjLf2hVuyKLxl1Jh/qqL4oEE4W5/Jfz64v9p6i+KOJ0CnMp0K/1xddui2PPwbz64nbVF0WcSmEuF9QjtjrJY+KpUSGSu2as4tXPdpCj+qKI4yjM5aLqVS3D/FFd6Ns6hlc/28ldb6zil5OqL4o4icJcfFI6wsvLN7fi/25sztK0AyQkprDpJ9UXRZxCYS4+M8YwpGNd5tzbiewci36TlzFn1Y92jyUiFEOYG2O6GWO+McZMMcZ0C/T9i/1a16nEorHxtK9XmUfe38Bj7+vNo0Xs5lOYG2OmG2MyjDGbzjve0xiz3RiTZox5LO+wBRwHIoH0wI4rTlGlbCneuKs9o7s3ZPaqH7l5ynLVF0Vs5Otf5jOAnvkPGGO8wESgFxALDDLGxALfWJbVC3gU+EvgRhWn8XoMD193Gf+8LY7vD56g94QUvlR9UcQWPoW5ZVlfA4fOO9weSLMsa7dlWVnAbKCPZVk5eV8/DJQK2KTiWNfk1Rd/Vz6SO2es4h+f7VR9UaSE+XPOPAbIf/UrHYgxxvQzxkwF3gImFHbjzMxM4uLifvtISkryYxSx22/1xSti+PtnO7hb9UWREhXmx21NAccsy7LmAfMuduPo6GhWr17tx8OL05SO8PLygFa0rluJp5M3k5CYwpQhbWkeU8Hu0URcz5+/zNOB2vk+rwXs9W8cCXbGGIaeX19crfqiSHHzJ8xXAY2NMfWNMRHAQGChrzc+cuQIw4cPJzk52Y8RxKl+rS+2q1eJR+Zu4PF5qi+KFCdjWRe/UGWMmQV0A6oC+4GnLMuaZoy5HngV8ALTLct6xtcHjouLs3Saxf2ycyxeXrydSV/uomWtCkwa3IZalaLsHkskKBljUi3Liivwa76EeXFQmIeWxZv38dCc9Xi9hvEDW9O1SbTdI4kEnQuFuV7OLyXi2st/x8KxufXF21//lvGfq74oEki2hbnOmYee+lXLMG9UZ268IoZXPt3BPW+u5sjJs3aPJeIKOs0iJc6yLN5asYe/LtpCjQqlmTykDZfXVH1R5GJ0mkUcxRjDbZ3q8e69ncg6l0O/Sct4T/VFEb8ozMU2bepUYtG4eNrWrcTDczfw+LyNnDmn+qLIpVCYi62qli3Fm3e1Z2S3hsz69gcGTFnOT7+csnsskaCjC6BiuzCvh0d7XsbUoW3ZnXmChPHf8PWOTLvHEgkqugAqjvLdgROMeCuVHRnHeLBHE0Z3b4THU9A2QCKhRxdAJWjUr1qG+aM706dVTV7+dAfDVF8U8YnCXBwnKiKMv99yBU/3uZyvd2bSe0IKm/fqzaNFLkRhLo70a31x9vB/1xffT9W7EIoURhdAxdHa1s2tL7apU4mH3lvPk/NVXxQpiC6ASlA4l53DS4u3M/Wr3bSqXZFJg9sQU7G03WOJlChdAJWgF+b18HivZkwZ0pZdGcdJGP8NKTsP2D2WiGMozCWo9Gz+OxaO6UJ0uVLcNn0lE5ekafdFERTmEoQaRJflX6O7kNCyJi99sp3hb63myCnVFyW0KcwlKEVFhPGPgVfw596xfLk9kxsmpLBl71G7xxKxjdosErSMMdzRpT7v3tuR02ez6Td5qeqLErLUZhFXyDx2hrGz1rBi9yGGdKzDnxJiKRXmtXsskYBSm0VcL7pcKWbe3YF7uzZg5oofGDB1BXu1+6KEEIW5uEaY18Pj1zdjypA2ufXFxBTVFyVkKMzFdXo2r8GCMV2oWjZC9UUJGQpzcaWG0WWZPyp/fTFV9UVxNYW5uFaZUrn1xad6x/Ll9gz6TEhh68+qL4o7qZoormaM4c4u9Zk9vCMns7LpO2kp89eqvijuo2qihIyMY6cZ+85aVn53iKEd6/KnhFgiwvR/TiV4qJooAlQrF8nb93RgeNcGvLViDwOmLufnI6ovijsozCWkhHk9PHF9MyYNbsPO/cdIGJ/CsjTVFyX4KcwlJF3fogYLxsRTqUwEQ6atZNKXadh1ylEkEBTmErIaVSvLgtFduL5FDV78OLe+ePS06osSnBTmEtLKlAojcVBr/jchliXbMrghMYVt+1RflOCjMJeQZ4zhrvj6zMqrL944cSn/WvuT3WOJFInCXCRPu3qVWTQunpa1KnL/u+v43wWbyDqXY/dYIj5RmIvk82t9cdiV9Xlz+R5uSVJ9UYKDXgEqcp5wr4cn/xDLpMFt2LEvr764S/VFcTa9AlTkAtIyjjNiZiq7M4/zSM/LuLdrA4wxdo8lIUqvABW5RL/WF3u1qMHzH21jxEzVF8WZFOYiF1GmVBgTBrXmf/7QjM+2ZtBnwlK27ztm91gi/0FhLuIDYwz3XNmAWcM6cvzMOW6cuJQF61RfFOdQmIsUQfv6lflgbDwtYipw3+x1/HnhZtUXxREU5iJFVK18JG8P68Dd8fWZsex7BiYtZ9+R03aPJSFOYS5yCcK9Hv6UEMuEW1uzbd8xEhK/UX1RbKUwF/FDQsuaLBzThQqlwxny2kqmfLVLuy+KLRTmIn5qVK0cC8bE06t5bn1x5Mw1HFN9UUqYwlwkAMqWCmPCrbn1xU+37qfPhKXs2K/6opQchblIgPxaX3znng4cO3OOPhNUX5SSozAXCbAODarwwdh4mseUV31RSow22hIpBtXKR/LOsI6/1RcH/XMF+4+qvijFRxttiRSzRRv28sjcDURF5J5X79igit0jSZDSRlsiNkpoWZMFo7tQvnQYg19bSdLXqi9K4CnMRUpA4+rlWDC6C9fGVufZD7cx6m3VFyWwFOYiJaRcZDiTBrfhyeubsXjLfvpMVH1RAkdhLlKCjDEM69qAt+/pwNFTubsvJq/fa/dY4gIKcxEbdGxQhQ/GxRNbozxjZ63lL8mbOZut+qJcOoW5iE2ql49k1vCO3NmlHq8v/Z6RM1MV6HLJFOYiNgr3eniq9+U83edyPtuawYNz1pOdo6aLFF2Y3QOICNzWqR4ns7J5/qNtlInw8ly/FnrjaCkShbmIQ4z4fUNOnjnH+C/SqBgVwWO9LrN7JAkiCnMRB3ngmiYcOpnFlK920TC6DDfH1bZ7JAkSOmcu4iDGGP7c+3KubFyVJ+ZvZOXug3aPJEFCYS7iMGFeDxNubUPtylGMmJnKnoMn7B5JgoDCXMSBKpQOZ/rt7bCAETPXcPpstt0jicMpzEUcql7VMrwyoBVbfz7K8x9ts3sccTiFuYiDXXVZ9d/2RF+8eZ/d44iDKcxFHO6Rnk1pHlOeR97fwN5fTtk9jjiUwlzE4UqFeUkc1Iaz53L443vrtRe6FEhhLhIE6lctwxN/aMayXQd5LzXd7nHEgYolzI0xZYwxqcaYhOK4f5FQNKhdHdrXq8wzH2zlwPEzdo8jDuNTmBtjphtjMowxm8473tMYs90Yk2aMeSzflx4F5gRyUJFQ5/EYnu3XnFNZ2TydvMXuccRhfP3LfAbQM/8BY4wXmAj0AmKBQcaYWGNMD2ALsD+Ac4oI0KhaOUZ1b8jC9XtZsj3D7nHEQXwKc8uyvgYOnXe4PZBmWdZuy7KygNlAH6A70BG4FRhmjCnwMTIzM4mLi/vtIykp6ZIXIRJKRnZrSKNqZXlqwWbOnNOLiSSXPxttxQA/5vs8HehgWdYYAGPMHcABy7IK3G0/Ojqa1atX+/HwIqGpVJiXp3rHMnTat7y5bA/DujaweyRxAH8ugBa02fJvnSnLsmZYlrXIj/sXkUJc2Tiabk2jGf/FTg6dyLJ7HHEAf8I8Hci/P2ctQO9MK1JCnry+GSezshn/+U67RxEH8CfMVwGNjTH1jTERwEBgoa83PnLkCMOHDyc5OdmPEURCV+Pq5RjUvjYzV+xhV+Zxu8cRm/laTZwFLAeaGmPSjTF3W5Z1DhgDfAJsBeZYlrXZ1weuUKECSUlJ9O7d+1LmFhHg/h5NiAz38tyH2ogr1Pl0AdSyrEGFHP8Q+DCgE4mIz6qWLcXIbg156ZPtpO45RNu6le0eSWyil/OLBLk7u9SjatlSvPDxdu3bEsJsC3OdMxcJjKiIMMZd3YhvvzvEVzsy7R5HbGLs+i95XFycpZ65SGBkncvh6le+pHxkOMlj4vF4CmoOS7AzxqRalhVX0Nd0mkXEBSLCPDx4TRM27z3Kh5t+tnscsYHCXMQlbmgVQ9Pq5Xh58Q7OZhf4wmtxMZ0zF3EJr8fwx+ua8t2BE8zVnuchx5+9Wfzya89cRAKnR7NqtKlTkX98tpO+rWOIDPfaPZKUEJ1mEXERYwwPX3cZ+46eZvrS7+weR0qQwlzEZTo1rEKPZtWZ+EUa+4+etnscKSEKcxEX+lNCM85mWzz/kV7mHyoU5iIuVLdKGYZ1rc/8tT/phUQhQm0WEZcae1VjGlUry2Pvb+Do6bN2jyPFzLYw166JIsUrMtzL325uxf6jp3l07gbt2+JyOs0i4mJX1K7IY70u46NN+5i4JM3ucaQY2dYzF5GSMezKBmzee5S/Ld5BtXKRDGhX++I3kqCjMBdxOWMML/ZvyaETWTw6bwMej6F/21p2jyUBptMsIiGgVJiXpKFxdG5YhT++t57xn+/kh4MndR7dRWzbArdx48ZW9+7d6d27ty6CipSQM+eyefDd9XywMXdnxT+0qMHfb7mCiDD9XRcMLrQFrvYzFwkxlmWxee9RPt60jwlL0rgmtjoTb22jQA8C2s9cRH5jjKF5TAX+eF1T/nLD5Xy6ZT+j3k7lzLlsu0cTPyjMRULY7Z3r8dcbm/PZ1gzGvLOW7BydQw9WCnOREDe0Y12e6h3Lp1v289TCTbooGqRUTRQR7uxSn31HTzP1q93UqFCa0d0b2T2SFJHCXEQAePS6y8g4eoaXPtlO9fKR6qIHGW20JSIAeDyGF25qSZdGVXhi3kbW/nDY7pGkCFRNFJH/8MvJLHpPSCHrXA7JY+OpVi7S7pEkj6qJIuKzilERJA2N4+ipc4yauYasczl2jyQ+UJiLyH9pVqM8L/Rvyeo9h/nroi12jyM+UJiLSIFuaFWTYVfW560Ve1iw7ie7x5GLUJiLSKEe6XkZcXUr8fi8jaRlHLd7HLkAhbmIFCrc6yHx1tZEhnsZ/fYaTmXpJf9OpTAXkQuqUaE0r95yBTsyjvH4vA3k6CX/jqQXDYnIRXVtEs1D1zThb4t3UDEqgru61Ld7pKDl9RpiKpYO+P0qzEXEJ6O7N+LwybNMS/mOGcu+t3ucoBVTsTRLH7sq4PerN6cQEZ9ZlsWS7RkcPnHW7lGCVlSEl14talzSbfXmFCIiLqBXgIqIuJzCXETEBRTmIiIuoDAXEXEBhbmIiAsEZZgnJSXZPUJAuGUdoLU4lVvW4pZ1QPGtxTFhnv8dhwr6d/5j/vyP4cs7GxX2PQUdP//YhT4/fy3+PqkluZaL/duu56SwrwXjWgL98wXBuZZAPycXmtOX73HSz9eFBGWYB+pxivo9xfHL5g8n/YD6w2lh7g/9fBV+3J8w95dbflcuxLYXDRljMoE9+Q5VAI5c4N/5j1UFDlziQ+e/n6J+T0HHzz92oc/PX4s/67jQnL58T1HXcrF/2/WcFPa1YFxLoH++IDjXEujn5EJz+vI9Tvr5qmtZVnRBX7AtzEVEJHAcc5pFREQuncJcRMQFFOYiIi7gqjA3xjQzxkwxxsw1xoy0ex5/GGNuNMb80xizwBhzrd3z+MMY08AYM80YM9fuWS6FMaaMMeaNvOdjsN3zXKpgfx7yc9nvR2Byy7IsR3wA04EMYNN5x3sC24E04DEf78sDTHPJWiq5aC1z7f45u5R1AUOB3nn/ftfu2f19fpz0PARgLbb+fgR4LX7llu2LzreQrkCb/IsHvMAuoAEQAawHYoEWwKLzPqrl3eYGYBlwa7CvJe92LwNtXLIWx4RIEdf1OHBF3ve8Y/fsl7oOJz4PAViLrb8fgVpLIHLLMW8bZ1nW18aYeucdbg+kWZa1G8AYMxvoY1nWc0BCIfezEFhojPkAeKcYRy5UINZijDHA88BHlmWtKeaRCxWo58VpirIuIB2oBazDYacmi7iOLSU8XpEUZS3GmK044PejMEV9XgKRW476wSxADPBjvs/T844VyBjTzRgz3hgzFfiwuIcroiKtBRgL9AD6G2NGFOdgl6Coz0sVY8wUoLUx5vHiHs4Pha1rHnCTMWYyULwv4wuMAtcRRM9DfoU9J07+/ShMYc9LQHLLMX+ZF8IUcKzQVzlZlvUl8GVxDeOnoq5lPDC++MbxS1HXchAIhl+4AtdlWdYJ4M6SHsYPha0jWJ6H/Apbi5N/PwpT2Fq+JAC55fS/zNOB2vk+rwXstWkWf2ktzueWdbllHaC1+MzpYb4KaGyMqW+MiQAGAgttnulSaS3O55Z1uWUdoLX4zu6rvvmu9M4CfgbOkvtfsLvzjl8P7CD3KvCTds+ptQTvWty4LresQ2vx/0MbbYmIuIDTT7OIiIgPFOYiIi6gMBcRcQGFuYiICyjMRURcQGEuIuICCnMRERdQmIuIuIDCXETEBf4faLjw7hKiS3kAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXMAAAD5CAYAAADV5tWYAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAbUUlEQVR4nO3deXRUVaL24d+uVAbGgISpAxiGyCSjQUZJGERQIzigTXNRFEyDiEiH1eLS1dpyb3s/DaDIICi0chFRQWVSUCEDYGSIFxCCDGGQADIJEQgQAuf7w+Cl6YRUqArnVOV91spa5FQq9W6Lejnu2rWPsSwLERHxby67A4iIiPdU5iIiAUBlLiISAFTmIiIBQGUuIhIAVOYiIgHAbdcDR0REWFFRUXY9vIiI38nIyDhmWVb1wm6zrcyjoqLYsGGDXQ8vIuJ3jDH7irpN0ywiIgFAZS4iEgBU5iIiAUBlLiISAFTmIiIBQGUuIhIA/K7ML17Slr0iIlfzqzLPy7/EfZNX89aKnZy7cNHuOCIijuFXZZ6bl8/N1coz/usd9JqYxopth+2OJCLiCH5V5lXKhzB14G3MGdKe4CDDkPc38MR769l77Izd0UREbOVXZX5Zl+gIvhzVlRfubsra3cfpNTGN15f/SG5evt3RRERs4ZdlDhDidvFk1wYkj4njnpa1mZKcRY/xqSzdfAhd11REyhrbyjwnJ4eEhAQWL17s1e+pUTmMiY+05uM/d6RK+RBGzP2ege+uZefhUz5KKiLifMaus9iYmBjL17sm5l+8xIfrfuL15dvJzbvI4E5RjOoZTaWwYJ8+joiIHYwxGZZlxRR2m99OsxTGHeRiUMcoksfE0T+mDjPX7KFbUioLMrK5pPXpIhLAAqrML6tWMZRXH2jJ5091JrJqORI/2UT/6elsOZBjdzQRkVIRkGV+Wau6VfhseCdee7Ale4+dIX7yal78/AdO5ubZHU1ExKcCuswBXC7Dw+3qsnJMHI91jGLu2p/olpTCB2v3aWsAEQkYAV/ml4WXC+bl+5qz9Jk7iK5ZiRc+20LfKavJ2HfC7mgiIl4rM2V+WdPalfkooQNv/rE1R0+d58Fp3zLmk00cPXXe7mgiItetzJU5gDGGvq0jWZEYx59jG7Bw4wG6J6Uwa/Ue8i9esjueiEiJlckyv6xiqJvn+zRl2bNdaV2vCq8syeSeSatJzzpudzQRkRIp02V+WcPqFZn9xO1MH3Qbp8/nM+Cd7xj54f9yKOes3dFERDyiMi9gjOGu5rVYkRjLqB7RLN/6Mz3GpzI1ZRfn87V3uog4m8r8KmHBQYy+8xa+GR1L50YRvLZsO73fWEXK9iN2RxMRKZLKvAj1qpXnnUdjeO/xdgAM/ud6npy9gf2/5NqcTETk36nMixHXuAbLnr2D53o3Yc2uY/SYkMqEr3dwNk9TLyLiHCpzD4S6gxge15AVibHc1bwWk1bspOeEVJZt+Vl7p4uII6jMS6B2eDneGtCGD5/sQMVQN8PmZPDorHVkHT1tdzQRKeNU5tehY8NqLH2mCy/FN2PjTyfp/UYar365jdPnddk6EbGHyvw6uYNcPN65PivHxNGvdSTTU3fTY3wKCzce0NSLiNxwKnMvVa8Uyuv9W/HpU52oUSmMUfM28siM79h26Fe7o4lIGaIy95G29ary+YjO/OP+Fuw4fIp7Jq3i5UVbyTl7we5oIlIGqMx9KMhl+FP7eiQnxvGn9vWYnb6X7kkpfLx+vy5bJyKlSmVeCqpWCOE/+7Vg0dNdiIqowF8XbOb+ad+yaf9Ju6OJSIBSmZeiWyPDmT+sIxMebsWBE2fpN3UNYxds5vhp7Z0uIr6lMi9lxhgeaFuH5DGxDOlcn/kZ2XRLSmF2+l7tnS4iPqMyv0EqhQXz4r3N+HLUHdwaGc7fFm4lfvIa1u/9xe5oIhIAVOY3WHTNSnwwtD1TB7YlJzeP/m+nM/qjjRz59Zzd0UTEj6nMbWCM4e4WtfkmMZanuzVi6eZDdEtKYUZaFnn5mnoRkZJTmduofIibMXc15qvRXWnfoBr/+OJH+ryZxuqdx+yOJiJ+RmXuAFERFZg1uB0zH4vhwkWL/5i5luFzMsg+ob3TRcQzbrsDyP/p0bQmnRtF8E7abqak7CJ5+xFGxDXiya4NCAsOsjueiDiYzswdJiw4iJE9olmRGEf3JjUY//UOek1MY8W2w3ZHExEHU5k7VGSVckwdeBtzhrQnOMgw5P0NPPHeevYeO2N3NBFxIJW5w3WJjuDLUV154e6mrNvzC70mpvH68h/JzdPe6SLyf1TmfiDE7eLJrg1YmRjLvS1rMyU5ix7jU1m6+ZD2ThcRQGXuV2pUDmPCI635ZFhHqpQPYcTc7xn47lp2Hj5ldzQRsZnK3A+1i7qJJSO7MK5vc7YcyKHPm6sYtySTX89p73SRssrnZW6McRlj/ssY85Yx5jFf/375TZDLMKhjFMlj4ugfU4dZa/bQPSmVBRnZ2jtdpAzyqMyNMbOMMUeMMVuuOt7bGLPdGLPLGDO24HBfIBK4AGT7Nq5crVrFUF59oCULR3SmTtVyJH6yif7T09lyIMfuaCJyA3l6Zv4e0PvKA8aYIGAK0AdoBgwwxjQDGgPplmX9BRjuu6hyLS3rVOHT4Z147aGW7D12hvjJq3nx8x84mZtndzQRuQE8KnPLstKAq/dqvR3YZVnWbsuy8oB5/HZWng2cKPiZi74KKsVzuQwPx9Rl5Zg4HusYxYfr9tMtKYUP1u7joqZeRAKaN3PmkcD+K77PLjj2KXCXMeYtIK2oOx89epSYmJjfv2bMmOFFFLlSeLlgXr6vOUtGdiG6ZiVe+GwLfaesJmPfieLvLCJ+yZu9WUwhxyzLsnKBIcXduXr16mzYsMGLh5fiNK1dmY8SOrBo00H+8cU2Hpz2LQ/dVofnejeheqVQu+OJiA95c2aeDdS94vs6wEHv4oivGWPo2zqSlYlxDIttyMKNB+ielMKs1Xt02TqRAOJNma8Hoo0x9Y0xIcAfgUW+iSW+ViHUzdg+TVj2bFda16vCK0syuWfSatKzjtsdTUR8wNOliR8C6UBjY0y2MWaIZVn5wNPAcmAb8LFlWVtLL6r4QsPqFZn9xO1MH3Qbp8/nM+Cd73h67vccyjlrdzQR8YKxa2+PmJgYS3Pm9jp34SLTUrJ4OzULlzGM7NGIIV3qE+rW3ukiTmSMybAsK6bQ2+wq8+joaKtbt27Ex8cTHx9vSwb5zf5fchm3JJOvMg9TP6ICL8U3I65xDbtjichVHFnmOjN3npTtR/j74kz2HDvDnc1q8rd7m1H3pvJ2xxKRAtcqc220Jb+La1yDZc/ewXO9m7Bm1zF6TEhlwtc7OJunz36JOJ3KXP5FqDuI4XENWZkYR+/mtZi0Yic9J6SybMvP2jtdxMFU5lKoWuFhTBrQhg+f7EDFUDfD5mTw6Kx1ZB09bXc0ESmEylyuqWPDaix9pgsvxTdj408n6f1GGq9+uY3T53XZOhEnUZlLsdxBLh7vXJ+VY+Lo1zqS6am76TE+hYUbD2jqRcQhbCvznJwcEhISWLx4sV0RpISqVwrl9f6t+PSpTtSoFMaoeRt5ZMZ3bDv0q93RRMo8LU2U63LxksXHG/bz2rIfyTl7gUc7RjH6zlsILxdsdzSRgKWlieJzQS7DgNvrkTwmjoHtb2Z2+l66J6Xw8fr9umydiA1U5uKVKuVDGNfvVhY93YWoiAr8dcFm7p/2LZv2n7Q7mkiZojIXn7g1Mpz5wzoy4eFWHDx5ln5T1zB2wWaOnz5vdzSRMkFlLj5jjOGBtnVYmRjL0C71mZ+RTbekFGan79Xe6SKlTGUuPlcpLJgX7mnGl6PuoEWdcP62cCvxk9ewfu/Vl5EVEV9RmUupia5ZiTlD2jN1YFtycvPo/3Y6oz/ayJFfz9kdTSTgaAtcuSFy8/KZmpzFjLTdBAcZRvWMZnCn+oS4dT4h4iltgSuOsffYGcYtyWTFj0doWL0Cf7/vVrpER9gdS8QvaJ25OEZURAVmDm7HzMdiyL9k8R8z1zJ8TgbZJ3Ltjibi19x2B5CyqUfTmnRuFMG7q3YzOXkXyduPMCKuEU92bUBYsC5bJ1JSOjMX24QFB/F092hWJMbRvUkNxn+9g14T0/gm87Dd0UT8jspcbBdZpRxTB97GB0PbE+J2MXT2Bh7/5zr2HDtjdzQRv6EyF8fo3CiCL0fdwYv3NGX93hPcNTGN15f/SG6e9k4XKY7KXBwlOMjF0DsasDIxlntb1mZKchY9xqeydPMh7Z0ucg0qc3GkGpXDmPBIa+YP60jV8iGMmPs9A99dy87Dp+yOJuJI+tCQON7FSxZz1/1E0vLtnDmfz2OdohjVM5rKYdo7XcoWfWhIAsIvZ/J4ffl25q3/iWoVQnm+TxPubxOJy2XsjiZyQ+hDQxIQbqoQwqsPtGDhiM7UqVqOxE820X96OlsO5NgdTcR2KnPxOy3rVOHT4Z147aGW7D12hvjJq3nx8x84mZtndzQR26jMxS+5XIaHY+qyckwcgztF8eG6/XRLSuGDtfu4qMvWSRmkMhe/Fl4umJfim7P0mS7cUrMSL3y2hb5TVpOx74Td0URuKJW5BIQmtSozL6EDkwa04eip8zw47VvGfLKJo6d02TopG1TmEjCMMdzX6g+sTIxjWGxDFm48QPekFGat3sMFXbZOApzKXAJOhVA3Y/s0YdmzXWlzc1VeWZLJvZNWk5513O5oIqVGZS4Bq2H1irz/eDtmDLqNM3n5DHjnO56e+z2Hcs7aHU3E51TmEtCMMfRqXotv/hLLsz2j+TrzMN2TUpmasovz+RftjifiM/o4v5Qp+3/JZdySTL7KPEz9iAq8FN+MuMY17I4l4hF9nF/kKqk7jvL3RVvZfewMdzaryd/ubUbdm8rbHUvkmvRxfpGrxN5SnWXPdmVsnyas2XWMHhNSmfD1Ds7maepF/JPKXMqsELeLYbENWZkYR+/mtZi0Yic9J6SybMvP2jtd/I7KXMq8WuFhTBrQhnkJHagY6mbYnAwenbWOrKOn7Y4m4jGVuUiBDg2qsfSZLrwc34yN+0/S+400Xv1iG6fP67J14nwqc5EruINcDO5cn+QxcdzfJpLpabvpnpTCwo0HNPUijqYyFylERMVQXnuoFZ891YmalcMYNW8jj8z4jm2HfrU7mkihVOYi19CmXlU+H9GZVx9owc7Dp7hn0ipeXrSVnLMX7I4m8i9U5iLFCHIZBtxej+QxcQxsfzOz0/fSPSmFj9fv55L2TheHUJmLeKhK+RDG9buVxSO7UD+iAn9dsJn7p33Lpv0n7Y4mojIXKanmfwjnk2EdmfhIKw6ePEu/qWsYu2Azx09r73Sxj8pc5DoYY7i/TR1WJsYytEt95mdk0y0phdnpe8nX3uliA220JeIDOw+f4uXFW1mz6zhNa1fmlb7NaRd1k92xJMBooy2RG8CyLJZt+ZlxSzI5mHOOfq3/wPN3N6Vm5TC7o0mA0EZbIjeAMYY+LWrzTWIsI7s34osffqZ7Ugoz0rLIy9fUi5QulbmIj5UPcZPYqzFfje5KhwbV+McXP9LnzTRW7TxqdzQJYCpzkVISFVGBmYPbMfOxGPIvWQyauY7hczLIPpFrdzQJQG67A4gEuh5Na9K5UQQzV+/hrZU7Sd5+hBFxjXiyawPCgoPsjicBQmfmIjdAWHAQI7o1YkViHD2a1GT81zvoNTGNbzIP2x1NAoTKXOQGiqxSjikD2/LB0PaEuF0Mnb2Bx/+5jj3HztgdTfycylzEBp0bRfDlqDt48Z6mrN97grsmpvH68h/JzdPe6XJ9VOYiNgkOcjH0jgasTIzl3la1mZKcRY/xqSzdfEh7p0uJqcxFbFajchgTHm7N/GEdqVo+hBFzv+dP76xlx+FTdkcTP6IyF3GImKibWDyyC+P63UrmoV/p8+Yqxi3J5Ndz2jtdiqcyF3GQIJdhUIebSR4Tx8MxdZm1Zg/dk1JZkJGtvdPlmlTmIg50U4UQXn2gBYtGdKFO1XIkfrKJ/tPT2XIgx+5o4lAqcxEHa1EnnE+Hd+L1h1qy7/gZ4iev5oXPfuDEmTy7o4nDqMxFHM7lMvSPqcuKxDgGd4pi3vr9dBufwgdr93FRUy9SQGUu4ifCywXzUnxzlj7ThcY1K/HCZ1voO2U1GftO2B1NHMC2Ms/JySEhIYHFixfbFUHELzWpVZl5CR14a0Abjp3K48Fp3zLmk00cPaXL1pVlujiFiB87cz6fycm7eHfVbsLcQYy+8xYGdbyZ4CD9T3cg0sUpRAJUhVA3z/VuwvJnu9Lm5qq8siSTeyetJj3ruN3R5AZTmYsEgAbVK/L+4+2YMeg2zuTlM+Cd73h67vccyjlrdzS5QVTmIgHCGEOv5rX45i+xjO55C19nHqZ7UipTU3ZxPv+i3fGklKnMRQJMWHAQo3pG881fYul6SwSvLdtO7zdWkbL9iN3RpBSpzEUCVN2byjN9UAzvP3E7Bhj8z/UMfX8DPx3XZesCkcpcJMDF3lKdZc92ZWyfJnybdYyeE1NJ3fHbxaU/WLuPZJ2xBwSVuUgZEOJ2MSy2ISsT4wBIKyjzaSlZDPufDLYe1J4v/k5lLlKG1AoPI+SqNejn8y8xbE4GOWe11a4/U5mLlHFNalXi0Mlz/G3hFrujiBdU5iJlXLM/VOaZHtEs3HiQRZsO2h1HrpPKXER4Kq4hretW4aWFW/hF2+v6JZW5iOAOcvH/HmzJqXP5/NfSbXbHkeugMhcRABrXqkRC1wYs+D6btbu1t4u/UZmLyO9Gdo+mdngY/7l0m6456mdU5iLyu3IhQfy1d2N+OJDDZ/97wO44UgIqcxH5F31bRdKqTjjjv9quDbr8iMpcRP6Fy2UYc1djDuac4+P1++2OIx5SmYvIv+nSKIJ2UVWZnLyLcxd0du4PVOYi8m+MMYy+8xYO/3qejzfo7NwfqMxFpFAdG1SjTb0qvLtqD/kXL9kdR4qhMheRQhlj+HPXhvz0Sy7Ltv5sdxwphspcRIp0Z7OaNIiowIy03ViW1p07mW1lnpOTQ0JCAosXL7YrgogUI8hleLxzFJuzc9i4/6TdceQabCvz8PBwZsyYQXx8vF0RRMQD97etQ8VQN/+Tvs/uKHINmmYRkWuqGOrmgbaRLNl8iGOnz9sdR4qgMheRYj3Sri55Fy+R/KOuF+pUKnMRKVatymEA5ObpA0ROpTIXkWKFuH+rirx8rTd3KpW5iBTr9zLXh4ccS2UuIsUKCfqtKs7rzNyxVOYiUixjDCFul7bEdTCVuYh4JDTIpTlzB1OZi4hHQtwqcydTmYuIR1TmzqYyFxGPhLpdegPUwVTmIuIRnZk7m8pcRDwS4nZpnbmDqcxFxCMhWs3iaCpzEfFIqDtI68wdTGUuIh7RnLmzqcxFxCMhWs3iaCpzEfGI3gB1NpW5iHgk1O3i/AWVuVOpzEXEI6E6M3c0lbmIeERLE51NZS4iHtFqFmdTmYuIR7TO3NlU5iLikRC3i0sW5Gve3JFU5iLiEV0H1NlU5iLikcvXAdW8uTOpzEXEI6HBuqizk6nMRcQjOjN3NpW5iHjk8py5zsydSWUuIh4JdevM3MlU5iLikVB3EIDWmjuUylxEPBKiM3NHU5mLiEe0ztzZVOYi4pHLc+baBteZVOYi4hGdmTubz8vcGBNnjFlljHnbGBPn698vIvbQOnNn86jMjTGzjDFHjDFbrjre2xiz3RizyxgztuCwBZwGwoBs38YVEbvoDVBn8/TM/D2g95UHjDFBwBSgD9AMGGCMaQassiyrD/Ac8HffRRURO2lporN5VOaWZaUBv1x1+HZgl2VZuy3LygPmAX0ty7r8z/YJINRnSUXEVvoEqLO5vbhvJLD/iu+zgfbGmAeAu4AqwOSi7nz06FFiYmJ+/z4hIYGEhAQv4ohIaQrVG6CO5k2Zm0KOWZZlfQp8Wtydq1evzoYNG7x4eBG5kfQGqLN5s5olG6h7xfd1gIPexRERp3K5DMFBRtMsDuVNma8Hoo0x9Y0xIcAfgUW+iSUiThQSpIs6O5WnSxM/BNKBxsaYbGPMEMuy8oGngeXANuBjy7K2ll5UEbFbiFtl7lQezZlbljWgiONfAF/4NJGIOJbK3Lls+zh/Tk4OCQkJLF682K4IIlJCoe4grTN3KG9Ws3glPDycGTNm2PXwInIdQtwuLU10KG20JSIe0xugzqUyFxGPhbhdWproUCpzEfFYqMrcsVTmIuIxrWZxLpW5iHgsVGXuWFqaKCIe02oW59LSRBHxmNaZO5emWUTEY1qa6FwqcxHxmN4AdS6VuYh4zGV+u8ivOI/KXEQkAKjMRUQCgMpcRCQAaJ25iEgA0DpzEZEAoGkWEZEA4FdlHkhn8hqLM2kszqSxFM/2Mr9yzry4P3v7H8GT+fmifqaw41cf01iuT1kay7Vud9JYvHlOirpNY/HN37Gi+FWZ+/KxSvozvn5CvaWxFH3cyWO5kaXhDacVoDcCaSzXYizLns9zGWOOAvuAcCCn4HBxf44AjnnxsFf+zpL+TGHHrz6msVyfsjSWa93upLF485wUdZvG4v1YbrYsq3phN9hW5iIi4ju2T7OIiIj3VOYiIgFAZS4iEgACpsyNMU2NMW8bY+YbY4bbnccbxph+xph3jDELjTG97M7jDWNMA2PMTGPMfLuzXA9jTAVjzPsFz8dAu/N4w9+fiysF2GvEN91lWZbtX8As4Aiw5arjvYHtwC5grIe/ywXMDJCxVA2gscy3++/Z9YwLGATEF/z5I7uz++I5ctJz4YOx2Poa8fFYvOou2wddMIiuQNsrBw4EAVlAAyAE2AQ0A1oAS676qlFwn/uAb4E/+ftYCu43HmgbIGNxTIGUcFzPA60Lfmau3dm9GYsTnwsfjMXW14ivxuKL7rJto60rWZaVZoyJuurw7cAuy7J2Axhj5gF9Lct6Fbi3iN+zCFhkjFkKzC3FyEXyxViMMQb4b+BLy7K+L+XIRfLV8+I0JRkXkA3UATbiwGnJEo4l8wbHK5GSjMUYsw0HvEaKUtLnxRfd5bi/nFeIBPZf8X12wbFCGWPijDGTjDHTgS9KO1wJlWgswEigJ/CQMWZYaQa7DiV9XqoZY94G2hhjni/tcF4oalyfAg8aY6YB/rJfc6Fj8aPn4kpFPS9Ofo0UpajnxSfd5Ygz8yKYQo4V+Qkny7JSgJTSCuOlko5lEjCp9OJ4paRjOQ74w4ut0HFZlnUGePxGh/FSUWPxl+fiSkWNxcmvkaIUNZYUfNBdTj4zzwbqXvF9HeCgTVm8pbE4XyCNS2NxplIdi5PLfD0QbYypb4wJAf4ILLI50/XSWJwvkMalsThT6Y7F7nd9C97J/RA4BFzgt3+9hhQcvxvYwW/vAL9gd06NxX/HEqjj0lic+WXHWLTRlohIAHDyNIuIiHhIZS4iEgBU5iIiAUBlLiISAFTmIiIBQGUuIhIAVOYiIgFAZS4iEgBU5iIiAeD/A8qBw9cVIRZaAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAA060lEQVR4nO3deVhV1f7H8fdiVBBRmVRQESHFGSQ1zQGHUhNtMNNs1CIr01u3wea61bW5xCEztbLJytvgPITzkPM8MjghCiiKKDKv3x/H7uVHoAgH9j6H7+t5eB7OBvb+Lg/n42KdtddSWmuEEELYFwejCxBCCGF9Eu5CCGGHJNyFEMIOSbgLIYQdknAXQgg7JOEuhBB2SMJdCCHskIS7EELYISdrn1ApFQS8DHhqrYdcOdYNGHHlei211l2sfV0hhBD/o8pyh6pSahYwEEjVWrcucrwfMBFwBGZord8t8rW5f4V7kWO3A35a68+LX8Pb21sHBgaWsxlCCFE9bdu27YzW2qf48bL23L8CJgOz/zqglHIEpgB9gSRgi1JqntZ6/1XOcy/wSElfCAwMZOvWrWUsRwghBIBS6lhJx8s05q61XgOkFzvcEYjXWidqrXOBOcDgqxTQGMjQWl8oW8lCCCHKqyJvqPoDJ4o8TgL8lVJeSqlpQJhS6sUiXx8FfFnaydLS0oiIiPjvx/Tp0ytQmhBCVG8VeUNVlXBMa63PAqNL+MLrVzuZi4sL4eHhREVFERUVVYGyhBBCVCTck4BGRR4HAMnlPZmnp6f01oUQwkoqMiyzBQhRSjVVSrkAw4B55T1ZRkYG0dHRzJ8/vwIlCSGEgDL23JVSPwA9AW+lVBLwutZ6plJqDLAUy1TIWVrrfeUtRHruQghhPWWa514VQkJCdGRkpIy5CyHEdVBKbdNaRxQ/bvU7VMtLeu5C2Ja4lEz+s/0k+5IzCKjrxvj+LfCs6Wx0WeIK6bkLIcpMa83GxLN8tiqBtXFncHZUBPt6cDglk/q1a/DpsPbcGFjP6DKrldJ67qYJ94iICC13qAphToWFmmX7U/hsdQK7TpzHu5YrD3Vpwr2dmlDP3YUdx88xbs5Oks5l8VSvEJ7qFYyTo6xLWBUk3IUQ162wUDN/dzIxsXEkpF2icT03HusRxF3hAdRwdvx/35uZncfr8/bxy/aTdGhSl0/vaU+jem4GVV59mD7cZVhGCPPQ2tJT/3jZYQ6lZNKivgdPRgbTv3X9a/bIf995kld+3QvAO3e2YVC7hlVRsu06mwD1gkCVdF/otZk+3KXnLoTxtNasjTvDR8sOsSspgyBvd57uewO3tWmAg0PZw+dEehbj5uxg+/Hz3Bnuz78Gt6aWq2nmb5hD3mVY/R5smAR3fA5thlz7Z0pg+tkyQghjHDt7idkbj5GVm09C6iU2H03Hv05N3h/SljvD/Ms1dt6onhs/PXYTMSvimbwijm3HzjFxWBjtG9WxfgNsUcJKWPA0nDsCYfdBs15Wv4Rpeu4yLCNE1UpMu8iUlQn8tvMkjg4Kz5rOuDo5EN09iHtubISrk+O1T1IGW46m8485O0m5kM3TfW9gdI9mOF7HXwF25dJZWPYK7PreMhQTNRGadq/QKWVYRggBQHxqJpNXxDNvVzIuTg6M6NSEx7oH4Vu7RqVdM+NyHi/9uoeFu0/ROagen9zTngaeNSvteqajNez+CZa+CNkZ0HUcdH8OnCv+byDhLkQ1dzglk0kr4lmwO5kaTo48cFMTHukWhI+Ha5VcX2vN3G1JvD5vH86ODrx3Vxv6tW5QJdc2VPoRWPgMJKwA/wgYFAN+rax2ehlzF6KaOnDqApNWxLFoz2ncXRwZ3aMZj9zcFK9aVRPqf1FKcXdEIyIC6zFuzg5Gf7udezs15rWBLf82rdIuFOTDn1Ng5QRwcIIBH0LESHComraapucuY+5CWNfekxnExMaxbH8KHq5OPNQ1kJFdm1LX3cXo0sjNL+SjZYf4fE0ioQ1qM3VEOE293Y0uy3pObof5Y+H0Hmg+wBLsnv6VcikZlhGimth14jyTVsTxx4FUPGo4MbJrU0Z2bYqnm/nWfVlxMIVnftpFXn4h797VlihbnxOfcxFW/hs2fQbuvjDgAwiNKvcc9rKQcBfCzm0/fo6Y2DhWHUrDs6Yzj9zclAe7BlK7hvlCvajk85d56ocdbDt2jhGdGvOqrQ7THF4GC/8JGcchYhT0eR1qeFb6ZWXMXQg7tfVoOhNj41gbd4a6bs48368593dugofJQ/0vDevUZE50Zz5cdojPVyey4/h5po4IJ9BWhmkupsKS8bD3P+DdHEYuhcadja5Keu5C2Ko/E88SExvHhoSzeLm7EN09iPs6N8Hdhu8E/WuYpqBQEzMsjMgWvkaXVDqtYcc3lnnreZeh27Nw8z/AqcrfqDb3sIy8oSrEtWmt2Zhwlk9j49h8JB0fD1ce6x7EvZ0a4+Ziu6FeVNK5LB77Zhv7T13g2Vua80TPZqhKHLMulzPxsOAfcHQtNO5iuRnJ5wZDSjF9uEvPXYjS/bXmS0xsHFuPncOvtiujezRjeMfGtjk+fQ2XcwsY/8tuft+ZTP/W9fnw7nbm+IskPxfWT4Q1H4BzDej7FoTdDw7GLW8sY+5C2CCtNasOpxETG8eO4+dp4FmDtwa34u6IRnYZ6n+p6eLIp/e0p3VDTyYsPkBi2iVmPBhh7BLCJzbDvLGQdgBa3QH93gMPP+PquQbpuQthQlprYg+kErMijt1JGfjXqckTkc0Y0iHAamu+2Ip1cWd48vvtODsqvngggrDGdau2gOwMiP0XbJkJtf3hto+geb+qreEqZFhGCBvw145Hk1bEsS/5Ao3q1WRMZDB3hAXg4lR9dzaKT73IyK+2kHIhm4nD2lfdsgUHFsCiZyHzNHQaDb1eAddaVXPtMpJwF8LECgs1S/adJiY2joOnMwn0cmNMrxAGt2+Is2xXB8DZizk8MnsrO0+c56X+oTzSrWnlvdF6IRkWPQcHF4Bfa4iKgYAOlXOtCqqyMXelVBDwMuCptR5y5ZgD8BZQG9iqtf7a2tcVwhYVFGoW7TnFpBVxHE65SJCPO5/c046otg1lD9JivGq58sOjnXnmp528s+gAJ85l8XpUK+suH1xYCNtmwR9vQkEu9HkDbhoDjrZxz0BRZQp3pdQsYCCQqrVuXeR4P2Ai4AjM0Fq/q7VOBEYppeYWOcVgwB9IB5KsVbwQtiq/oJAFuy2hnpB2iRDfWsQMD+O2Ng2q71rnZVDD2ZHJw8N5t+5Bpq9J5HxWHh8NbWedv25SD8D8cXBiEwT1hIGfWNZct1Fl7bl/BUwGZv91QCnlCEwB+mIJ7C1KqXla6/0l/HxzYKPW+vMroR9boaqFsFH5BYX8tjOZKSvjOXLmEi3qezDl3nD6t65/XdvYVWcODoqXBoRS182F95YcJDM7j6kjOlDTpZxvNOdlw9oPYd2n4Oph2fKu7T2Vuh5MVShTuGut1yilAosd7gjEX+mpo5Sag6WHXlK4JwG5Vz4vKF+pQtiuvIJCft1+kskr4zmenkXLBrWZdl8HbmnpJ6FeTo/3bIZnTWde/m0PD87azIyHIq5/HZ2j6yy99bPx0HYY3PoOuHtXTsFVrCJj7v7AiSKPk4BOSikv4B0gTCn1otZ6AvALMEkp1Q1YU9LJ0tLSiIj433sC0dHRREdHV6A8IYyXm1/I3G1JTF0VT9K5y7Tx9+SLByLoE+prvrsubdC9nRrjUcOJp3/cyfDpfzJ7ZMeyrVOflQ7LX7MsH1A3EO7/tVL2MTVSRcK9pN9MrbU+C4wudjALGHW1k/n4+CCzZYS9yMkv4KetSXy2Mp7kjGzaNarDW4Nb07O5j4S6lUW1a0itGk48/u02RszYxPePdqZeaWvWaw37foHFL1gCvus46DEeXAy8OaqSVORdiCSgUZHHAUByeU+WkZFBdHQ08+fPr0BJQhgrO6+Ar9Yfocf7q3j1t700qFOT2SM78tsTXYhsIb31yhLZ3JeZD97IkTOXGDFjE+cu5f79m84fh++HwtyR4BkA0aug77/sMtihYj33LUCIUqopcBIYBtxrlaqEsDGXcwv4btMxPl+TSFpmDh0D6/HR0HZ0aeYlgV5FugZ788UDETwyeyv3zdzEd490oo6bCxQWwKbPYcXblm+8dQJ0eqzKtrszSpluYlJK/QD0BLyBFOB1rfVMpdQA4FMsUyFnaa3fKW8hchOTsEVZufl8++cxpq9J5MzFXG4K8mJcnxA6B3kZXVq1tepQKtGzt9G8vgffR7nhsewZSN4BIbdYlg6o09joEq3K9HeoypK/wpZczMln9sajzFh7hPRLuXQL8eapXiF0bFrP6NIEsGrvUQ7/+AqjHBei3L1x6P+eZbEvO/wryvSrQnp6ejJ9+nSjyxDiqi5k5zF7w1FmrDvC+aw8ejb34aleIXRoUsWLWYnSxcfS84+n6el4jB8LIllRZwyTWvTCxQ6D/WpME+5/vaEqPXdhRhmX8/hy/RFmrTvChex8erfwZWzvENo1qmN0aeIvl87A0pdg94/gFQwPLUSdacLSubt5fu4uPh7avlrdU2CacJeeuzCj81m5zFp3hC/XHyUzJ59bWvoxtncIrf0rf+NjUUZaw645lmDPyYTuz0O3f4JzDYYGQlpmDh8sPYRv7Rq8NCDU6GqrjGnCXQgzSb+Uy4y1iXy94SiXcgvo37o+Y3oF06qhhLqppCfC/H/AkdXQqJNluzvf/x/gT/RsRsqFbKavScTXw5VHutnuejHXwzThLsMywgzOXMzhi7WJfLPxGJfzCritTQOe6hVC8/oeRpcmiirIg42TYdW74OhimQXTYWSJ290ppXg9qhVnLubw9sID+Hi4Mri9vwFFVy3ThLsMywgjpWZmM311It9uOkZufiGD2jVkTK9ggn0l1E0naRvMHwspeyE0Cvq/D7UbXvVHHB0UHw9tz5mLm3lu7m4a13Or+h2dqphppkLKPHdhhJQL2Xy2KoEfNh8nv1AzuH1DnowMppmPuXbbEVjG01e8A5umgUcDGPABhA68rlOcu5TL4CnrycotYN6YrjSsU7OSiq06Ms9diCKSz19m2uoE5mw5QUGh5q5wf57oGUygt7vRpYmSHFoCC/8JF07CjY9A79egRu1ynepwSiZ3Tt1AEy83fh59E24uphnAKBfTh7v03EVVOJGexWerE/h5q2VB0yEdAniiZzCN6tnn+iI2L/O0ZZGv/b+BTygMioFGHSt82pUHUxn59RYGtGnA5OFhNr1EhOlvYhKiMh0/m8WUlfH8Z3sSDkpxz42NGN2jGQF1JdRNqbAQdsyGZa9BfrZlY+ou48CplNUer1NkC1+ev7UF7y05SHjjuoy6ualVzmsmEu7Crh05c4kpK+P5dcdJHB0UIzo1ZnTPZjTwtP2xVruVdhgW/AOOrYfAbjDwU/AOtvplRvcIYsfxc0xYdIC2AZ7cGGhfS0eYZlhGxtyFNSWkXWTKinh+23kSZ0cHRnRqwmM9gvCrXcPo0kRp8nMsW92t/RCc3eCWtyHsvkpdD+ZCdh6DJq0jK7eAhWO74eNRho0+TEbG3EW1EJeSyaQV8czfnUwNJ0fu69yYR7sH4eshoW5qxzZatrs7cwhaD4F+E6CWb5Vc+sCpC9wxdT3hjevyzahONrdBuYy5C7t28PQFJsXGs2jvKWo6O/JY92Y80q0p3mXZck0Y5/J5iH0Tts4Cz8YwYi6E9K3SEkIb1ObNQa144T97mL4mkcd7NqvS61cWCXdh0/YlZzApNp4l+05Ty9WJJ3sGM/LmpqVvsybMQWs4MA8WPQ+XUuGmMdDzRXA15v6CoRGNWHP4DB8tO0SXZl52sSCchLuwSXuSMpgYG8cfB1LwqOHE2N4hjOwaaNl5R5hbRhIseg4OLYL6bWD4D+AfbmhJSin+fUcbdhw/x7g5O1g4thvurrYdj7Zdvah2dhw/x6QV8aw4mIpnTWee6XsDD3YJxLOms9GliWspLIAtMy3DMIUF0Pct6PwEOJojhjzdnPl0WBjDpm/kjXn7+ODudkaXVCHm+FdFFg4TV7ftWDoTY+NZcziNOm7OPHdrcx64qQkeNSTUbULKPssbpklboFkvGPgJ1A00uqq/6di0Ho/3bMaUlQkMaNOAyBZV86ZuZZDZMsLUNh9JZ2LsYdbHn8XL3YVHuwdxX+cm1LLxP5mrjbzLsOYDWD8RanhCv3ehzd2m3u4uJ7+AqEnruHA5n6VPdzf9X4UyW0bYDK01GxPPEhMbx5+J6XjXcuXlAaGM6NzY5tcBqVYSV1tuRkpPhPYjLPPW3cx/o5CrkyMf3t2OO6Zu4J2F+3l/iG0Oz8grRZiG1pr18ZZQ33w0HV8PV14b2JLhHRtT08XR6PJEWWWlw7JXYee3ULcpPPA7BPU0uqrr0jagDo91D2LqKsvwTM/mtjc8I+EuDKe1ZvXhNGJi49h+/DwNPGvwr8GtGBrRiBrOEuo2Q2vYMxeWjIfs83DzM9DjeXC2zaUexvUJYfn+FF76ZQ/Ln+lhc7NnrF6tUioIeBnw1FoPuXKsJ/AWsA+Yo7VeZe3rCtujtWbFwVRiYuPYlZSBf52avH17a+6OCMDVSULdppw7CguegYRY8O8AUb9D/dZGV1Uhrk6OTLizDUOmbSQmNo4XbWz/1TKFu1JqFjAQSNVaty5yvB8wEXAEZmit39VaJwKjlFJzi5xCAxeBGkCStYoXtklrzfL9KcSsiGPvyQsE1K3Ju3e24c7wAFyc/r5NmjCxgnzY9Bms/DcoB8uuSDc+Ag728Z9zRGA9hkYEMHPdEe4MD7Cp7RbL2nP/CpgMzP7rgFLKEZgC9MUS2FuUUvO01vtL+Pm1WuvVSik/4GNgRIWqFjapsFCzdN9pYlbEc+DUBZp4ufH+kLbcEeaPs6OEus1J3mnZ7u7ULrihP9z2IXgGGF2V1Y3vH8qy/Sm8+ttefnyss82s/V6mcNdar1FKBRY73BGIv9JTRyk1BxgM/C3ctdaFVz49B8hiH9VMQaFm8d5TTIqN51BKJkHe7nw8tB2D2jXESULd9uResvTU/5wK7j5w99fQcrCppzdWRD13F8b3a8H4X/bwy/aT3NXBNv4Dq8iYuz9wosjjJKCTUsoLeAcIU0q9qLWeoJS6E7gVqIPlL4C/SUtLIyLif1M1o6OjiY6OrkB5wmgFhZoFu5OZtCKe+NSLBPvWYuKw9gxs29DmVt4TV8T9AQuehozj0OFh6PMG1KxjdFWVbmhEI37aeoJ/LzrALa38bOLmuYqEe0mvTq21PguMLnbwF+CXq53Mx8cHuYnJPuQXFDJvVzKTV8STeOYSzf08mHxvGP1bN5BQt1UX02Dpi7DnZ/C+AR5eDE26GF1VlXFwULwxqBWDJq9n6qoEXujXwuiSrqki4Z4ENCryOABILu/JZPkB25dXUMivO04yZWU8x85m0aK+B5+NCOfWVvVxkFC3TVrDzu9g6cuW4Zge46HbM+BU/UZX2wbU4Y4wf2auO8KITo1Nv0VjRcJ9CxCilGoKnASGAfdapSphU3LzC/llexJTVsVzIv0yrf1rM/3+DvQJ9ZNQt2VnEyzrwRxdC41vgqiJ4NPc6KoM9dytzVm05xQfLD3ExGFhRpdzVWWdCvkD0BPwVkolAa9rrWcqpcYAS7FMhZyltd5X3kI8PT2ZPn16eX9cGCAnv4Cftybx2aoETp6/TLsAT94c1IrI5r42M6NAlKAgz7IWzOr3wamGZQ/T8AfBQd78blinJo92C2Lyynge7tqU9iZe9900C4fJHqq2IzuvgB+3nOCzVQmcvpBNWOM6jOsdQo8bfCTUbd2JLZbpjan7LTNg+r8PHvWNrspULubk0/ODVTT1duOnx24y/Hfe9AuHSc/d/LLzCvh+03GmrU4gNTOHGwPr8uHd7ega7GX4L7iooOwLsOIt2PwF1G4Iw36AFgOMrsqUark68Y8+Ibzy215WHU4j0qTrzkjPXVxTVm4+3/15nM/XJHLmYg6dg+oxtncINwVJqNuFgwth4bOQeQo6PQa9XgFX27kT0wi5+YX0+mgV9dxd+P3Jroa+DqTnLq7bpZx8Zm88xhdrE0m/lMvNwd481SuMTkFeRpcmrOHCKVj8vGUvU99WcM83EPC3jBAlcHFyYGzvEJ6fu5s/DqTSt6Wf0SX9jWnCXZhHZnYeszceY8baRM5l5dH9Bh/G9Q6mQxPzr8UtyqCwELZ9CX+8Afk50Ps16DIWHM1/Y46Z3Bnmz9SV8Xy8/DC9W/iabmaYacJd5rkbL+NyHl9vOMrMdUfIuJxHrxa+PNUrmLDGdY0uTVhL6kHL9MYTf0LT7paZMF7NjK7KJjk5OjCuTwhP/7iLpftO079NA6NL+n9MM+Yu2+wZ53xWLrPWH+XL9UfIzM6nb0s/xvYKoU2Ap9GlCWvJz4G1H8Haj8G1Ftz6b2g33G7Xg6kqBYWaWz5ZjaODYsm47ob03k0/5i6qXvqlXGauS+TrDce4mJNPv1b1eap3MK0aSqjblaPrLb31s3HQZqgl2Gv5GF2VXXB0UIztHcK4OTtZfiCFW1uZZ9qoaXruMlum6py9mMMXa48we+NRLucVMKBNA57qFUyL+rWNLk1Y0+VzsPx12P411GkMAz+B4D5GV2V38gsKifxoFd61XPnl8S5VPnPG9D13mS1T+dIyc5i+JoFv/zxOdn4BUW0bMqZXMDf4ybQ3u6I17PsVFr8AWWegy1PQ80VwcTe6Mrvk5OhAdLcgXv19H1uOnqNjU3NMPDBNuIvKk3Ihm89XJ/LdpmPkFRRye3t/nogMJti3ltGlCWs7fwIWPQuHl0CDdjDiZ2jY3uiq7N6QDo345I84pq1OkHAXle9UxmWmrUrghy0nKCjU3Bnmz5ORwQR6Sw/O7hQWwObpEPsWoOGWd6DTaHCUl3hVqOniyENdAvl4+WEOnc40xXZ8MuZuh06ev8zUlfH8vDWJQq0Z0iGAJ3oG09jL3EuUinI6vQfmjYXk7ZYx9ds+hrpNjK6q2jl3KZcu766gf5v6fDy0fZVdV8bcq4ET6VlMXRXP3G2WPcjvjmjEEz2bmX7daVFOeZdh1buwYRK41YO7ZkLru2R6o0HqurswrGMjvtl4jOdvbUF9zxqG1mOacBfld+zsJaasjOeX7SdxUIrhHRszukczGtapaXRporIkrLRsd3fuCITdB33fsgS8MNTDXZry1YajfL/5OM/0vcHQWiTcbVhi2kUmr4zn953JODko7r+pCY91b2Z4j0FUoktnYdnLsOsHqNcMHpxvudNUmEJjLzcim/vy/abjjIkMxsXJuDXwJdxtUHxqJpNWxDN/VzIuTg483CWQ6B5B+HpIqNstrWH3T5Z9TLMzoNuz0P05cJbn3Gzuv6kJD3+5hSX7TjOoXUPD6pBwtyGHTmcyaUUcC/ecoqazI492D+LRbkF416p++1lWK+lHLEMwiSsh4EaIigG/lkZXJUrRI8SHJl5uzN5wVMIdZOGwq9mffIFJK+JYvPc07i6OPN6jGY90C6Keu4vRpYnKVJAPGydb3jR1cIIBH0LESHBwNLoycRUODor7Ozfh7YUH2J98gZYNjbnz2zRTIWXhsL/bezKDmNg4lu1PwcPViYe7BjLy5qbUcZNQt3snt1u2uzu9B5rfBgM+AE9/o6sSZZSRlUenCX9wR5g/E+5sW6nXMv1USPE/O0+cZ1JsHLEHU6ldw4mn+9zAQ10D8awp623bvZyLsPId2DQN3H3hnm8hVP6StTWebs4MbufPrztOMr5/qCGvXQl3E9l27BwxsXGsPpxGHTdnnr3lBh7oEkjtGhLq1cLhZbDwGcg4ARGjoM/rUENW6LRVIzo35setJ5i/K5n7Olf9TWUS7iaw5Wg6MbFxrI07Qz13F17o14L7b2pCLVd5eqqFzBRYMh72/QI+LWDkUmjc2eiqRAW18fekuZ8HP29Lso9wV0oFAS8DnlrrIUWOuwNrgNe11gusfV1btDHhLDGxcWxMPIt3LRdeGtCCEZ2a4C6hXj1oDTu+gWWvWO42jXwZuo4DJ5n9ZA+UUtwdEcDbCw8Ql5JJSBWvvlqmGfZKqVlKqVSl1N5ix/sppQ4ppeKVUuMBtNaJWutRJZzmBeCnipds27TWrI8/w9DPNzL8iz+JT7vIqwNbsvb5XkR3bybBXl2ciYOvBsK8p8CvNTy+AXo8L8FuZ24P88fJQf13SZCqVNYk+QqYDMz+64BSyhGYAvQFkoAtSql5Wuv9xX9YKdUH2A9U2zsutNasiTtDTGwc246dw6+2K29EtWRYx8bUcJapbdVGfi6s/xTWfADONS1z1sPuBwfj7mQUlce7liuRLXz5ZcdJnru1OU6OVfc8lynctdZrlFKBxQ53BOK11okASqk5wGAsIV5cJOAOtAQuK6UWaa0Ly121DdFas+pQGhNj49h54jwNPWvw1u2tubtDgIR6dXN8k2V6Y9pBaHUn9HsXPPyMrkpUsrs7BLB8fwqrD6fRO7Tqnu+KjAH4AyeKPE4COimlvIB3gDCl1Ita6wla65cBlFIPAWdKCva0tDQiIv43VTM6Opro6OgKlGcsrTV/HEglJjaOPSczCKhbkwl3tuGu8ABD15sQBsjOgD/ehK2zwDMA7v0JbrjV6KpEFYls4YuXuwtztyXZTLiXtK6o1lqfBUaX9ANa669KO5mPjw/2cBNTYaFm2f4UYmLj2H/qAo3rufH+XW25I9wf5yr8k0yYxIH5sOg5uJgCnR+3vGnqKjtgVSfOjg7cHubP7I1HSb+UW2V3llckbZKARkUeBwDJ5T3ZX8sPzJ8/vwIlGaewULNw9ykGxKxl9LfbuJxXwEd3t2PFP3sw9MZGEuzVzYVkmDMCfrwP3LzhkT+g3wQJ9mrqrvAA8go0i/eeqrJrVqTnvgUIUUo1BU4Cw4B7y3syW92so6BQs3DPKSbFxhGXepFmPu58ek97BrZtUKVvngiTKCyErTMtwzCFedDnTbjpSXCUG9Gqs9AGHgT5uLNg1ylGdKqaOe9lCnel1A9AT8BbKZWEZa76TKXUGGAp4AjM0lrvK28htrZwWH5BIfN3JzNpRTyJaZe4wa8Wk4aHMaBNAxwdZCecaillP8wfB0mbIagnDPwE6gUZXZUwAaUUA9s2ZPKKOFIzs6tkee6yzpYZXsrxRcAiq1ZkcnkFhfy24yRTVsZz9GwWLep78NmIcG5tVR8HCfXqKS8b1n4I6z4FVw+443Noe49sdyf+n6i2DYiJjWPxntM82CWw0q9nmjtmzD4sk5tfyK87kpiyMoHj6Vm0alibz+/vQN9QPwn16uzIWktvPT0B2g2HW94Bdy+jqxImFOLnQXM/DxbsTq5e4W5WOfkFzN2WxNSVCZw8f5m2AZ68NjCC3qG+KOmZVV9Z6bD8NcvyAXUD4f7foFmk0VUJkxvYtgEfLT/M6YzsSt8O0zTv+Jlttkx2XgHfbDxK5AerePnXvfh4uPLlwzfy+5Nd6dPST4K9utIa9syFKR1h5/fQ9R/w+EYJdlEmA6/szLRwT+XPmjFNz90swzLZeQX8sPk401YnkHIhh4gmdXlvSFtuDvaWQK/uzh2Dhf+E+OXQMAzu+wUaVO5GDMK+NPV2p1XD2izYncyom5tW6rVME+5Gz5a5nFvAd5uOMW11Imcu5tCpaT0+Gdqem5p5SahXdwX5sPlzWPE2oCzLBnSMlu3uRLnc1rYB7y85RNK5LALqulXadUwT7kb13C/l5PPtn8f4Ym0iZy7m0qWZF5PvDaNzkLwpJoBTu2DeWDi1E0Juhds+gjqNrvljQpTmtjaWcF+2L4WRldh7N024V7WLOfnM3niUGWuPkH4pl24h3oztHcKNgfWMLk2YQW4WrJoAG6eAmxcM+RJa3SHTG0WFNfFy5wa/WsQerCbhXlXDMhey8/h6/VFmrj/C+aw8Ipv78FTvEMIb1620awobEx8LC56G88cg/EHo+ybUlN8PYT19Qv2YviaRjMt5lba/qmnCvbKHZTKy8pi1/giz1h8hMzufPqG+jO0dQtuAOpV2TWFjLp2BpS/B7h/BKwQeWgSBXY2uStih3qF+TF2VwOrDaQy6MoPG2kwT7pXl3KVcZq0/wlfrj5KZk8+trfx4qlcIrf1l42Fxhdaw6wdLsOdchB4vwM3PgHO13VtGVLL2jergXcuFP/anSLhfr7MXc5ix7gizNxwlK6+AAa0bMKZXMKENahtdmjCTswmw4B9wZA006gRRE8E31OiqhJ1zdFD0auHL4r2nySsorJRVY00T7tYac0/LzGHG2kS++fMYl/MKGNi2IU/1CuaGKt6cVphcQR5smASr3wNHF7jtY+jwsGx3J6pMn1A/ftqaxJYj6XQJ9rb6+U0T7hUdc0+9kM3naxL5btMxcvMLGdzenycjgwn2lfWzRTFJ2yzb3aXshdAo6P8B1G5gdFWimrk5xBtXJwf+OJBq3+FeXqczspm2OoHvNx+noFBze3t/noxsRpCPhLooJifTciPSps/BowHc8x2EDjS6KlFNubk40TXYm+UHTvPqwFCr3yxp8+H+2ap4vtt0nLvCA3gishlNvNyNLkmY0aHFlqUDLiRDx0eh16tQQ95/EcbqE+rHioOpxKVetPrQsc2H+5O9gnmkWxCN6lXebbzChmWehsUvwP7fwCcURn0FjToaXZUQAPQO9aXuUmeOnrkk4V5cVexoImxQYSFs/xqWvw752dDrFegyDpyqZnNiIcrCr3YNtr7St1J2bzPN1ACzLfkrbFjaIfhqgGWKY4O28PgG6P6cBLswpcraltM0PXezLPkrbFh+Dqz7BNZ+BM5uMHgKtB8h68GIask04S5EhRzbaNnu7swhaD3EsixvLR+jqxLCMBLuwrZdPg9/vAHbvgTPxjBiLoT0NboqIQwn4S5sk9aw/3dY/DxcSoObxkDkS+AiU2GFAAl3YYsykmDRc3BoEdRvC/f+aNn2TgjxX1YPd6VUEPAy4Km1HnLlWCgwDvAGYrXWn1n7uqIaKCyALTMg9l+Wz/u+BZ2fAEfpowhRXJmmQiqlZimlUpVSe4sd76eUOqSUildKjQfQWidqrUcV/T6t9QGt9WhgKBBhreJFNXJ6L8zsaxmGadQJnvwTuo6VYBeiFGWd5/4V0K/oAaWUIzAF6A+0BIYrpVqWdgKl1CBgHRBbrkpF9ZR3Gf54E6b3gHPH4M4ZcN9/oG6g0ZUJYWpl6vZordcopQKLHe4IxGutEwGUUnOAwcD+Us4xD5inlFoIfF/uikX1kbjaciNSeqJlvvotb4Ob7HErRFlU5G9af+BEkcdJQCellBfwDhCmlHpRaz1BKdUTuBNwBRaVdLK0tDQiIv43YhMdHU10dHQFyhM2Kysdlr0CO7+DekHwwDwI6mF0VULYlIqEe0m3/Wmt9VlgdLGDq4BVVzuZi4sL4eHhlb5BtjAxrWHPz7BkPGRnWLa66/E8ONc0ujIhbE5Fwj0JaFTkcQCQXN6TyfID1dy5o7DgGUiIBf8OEBUD9VsbXZUQNqsiC4dtAUKUUk2VUi7AMGBeeU8mC4dVUwX5sD4GpnSGE5ug//swarkEuxAVVKaeu1LqB6An4K2USgJe11rPVEqNAZYCjsAsrfW+8hYiPfdqKHkHzBsLp3fDDf3htg/BM8DoqoSwC0prbXQNAISEhOjIyEgZc68Oci/Byn/Dn1PB3QcGfAChg2T1RiHKQSm1TWv9t/uHTHMHiPTcq4m45Zax9Yzj0OFh6PMG1KxjdFVC2B3ZrENUjYupMHcUfDfEMvvl4SUQ9akEuxCVRHruonJpDTu+tcxbz8uCni/CzU+Dk6vRlQlh16TnLirPmXj4OgrmjQHfUBi9DnqOl2AXogpIz11YX34ubJgIqz8Apxow8FMIfxAcTNOXEMLumSbchZ04sQXmj4XU/dDyduj/HnjUN7oqIaod03SlZFjGxmVfgIXPWpblzc6A4XNg6NcS7EIYxDQ9dxmWsWEHF1qCPfMUdHoMer0Crh5GVyVEtWaacBc26MIpWPwcHJgPvq3gnm8hoIPRVQkhkHAX5VFYCNtmWTbRKMiF3q9Dl6fA0dnoyoQQV8iYu7g+qQfgy36w8J+WTakf3wDdnpFgF8JkTNNzlzF3k8vLhrUfwbpPwLUW3P4ZtBsu68EIYVKmCXdhYkfXw/xxcDYO2t4Dt/4b3L2NrkoIcRUS7qJ0l8/B8tdg+2yo0wTu+wWCextdlRCiDCTcxd9pDft+hcUvQNZZ6DLWsmyAi7vRlQkhykjCXfx/509Y3iyNWwoN2sN9c6FBO6OrEkJcJ5ktIywKC+DPz2BKJzi61jKu/kisBLsQNso0PXeZLWOgU7st68Ek74DgvnDbR1C3idFVCSEqwDThLgyQmwWr34UNk8GtHtw1E1rfJdMbhbADEu7VVcIKWPA0nDsKYfdD339ZAl4IYRck3KubS2dh6Uuwew54BcODC6BpN6OrEkJYmYR7daE17P4RlrwIOReg+3PQ7VlwrmF0ZUKISlAp4a6UCgJeBjy11kOuHLsduA3wBaZorZdVxrVFCdITLUMwiasg4EaIigG/lkZXJYSoRGWeCqmUmqWUSlVK7S12vJ9S6pBSKl4pNR5Aa52otR5V9Pu01r9prR8FHgLusULt4loK8mDdpzC1CyRtgwEfwshlEuxCVAPX03P/CpgMzP7rgFLKEZgC9AWSgC1KqXla6/1XOc8rV35GVKaT22DeOEjZAy0GwoAPoHZDo6sSQlSRMoe71nqNUiqw2OGOQLzWOhFAKTUHGAz8LdyVUgp4F1istd5e7orF1eVchBVvw+bPoZafZQON0CijqxJCVLGKjrn7AyeKPE4COimlvIB3gDCl1Ita6wnAU0AfwFMpFay1nlb0RGlpaURERPz3cXR0NNHR0RUsr5o5vNSydEBGEtw4Cnq/BjU8ja5KCGGAioZ7SXe7aK31WWB0sYMxQExpJ3JxcSE8PJyoqCiioqSneV0yU2DJC5bFvnxCYeRSaNzJ6KqEEAaqaLgnAY2KPA4AkstzIll+oBwKC2HHN7D8Vci7DJGvQNdx4ORidGVCCINVdOGwLUCIUqqpUsoFGAbMK8+JZOGw65R2GL4eaFkTxq+1Zbu7Hs9JsAshgOvouSulfgB6At5KqSTgda31TKXUGGAp4AjM0lrvK08h0nMvo/xcWP8prPkAnGvCoEnQ/j5wMM0Cn0IIE1Baa6NrACAkJERHRkbKmPvVHP/Tst1d2kFodSf0exc8/IyuSghhIKXUNq11RPHjpll+QHruV5GdAX+8CVtngmcjuPcnuOFWo6sSQpiYaf6WlzH3UuyfB5M7wrYvofOT8MSfEuxCiGuSnrtZZZyExc/DwQXg1waGfw/+HYyuSghhI6TnbjaFBbD5C8t2d/GxlnXWo1dKsAshrov03M0kZZ/lDdOkLRAUCQM/hnpBRlclhLBBpgn3ai0v2zK1cf2nluUC7pgObYfKdndCiHKTcDfakTUw/x+QngDthsMt74C7l9FVCSFsnIy5GyUrHX5/Er6OAl0A9/8Gd0yTYBdCWIVpeu7VZsxda9j7H1gy3hLwNz8N3Z8HFzejKxNC2BHThHu1cO6YZUne+OXQMBzu/xXqtzG6KiGEHZJwrwoF+bBpGqx8B1CWZQM6RoODo9GVCSHslIy5V7bknTCjFyx7GQK7wZOboPPjEuxCiEplmp673Y25516CVRNg41Rw84K7v4KWt8v0RiFElTBNuNuV+D9gwdNw/jiEPwh934SadY2uSghRjUi4W9PFNFj6Euz5CbxC4KFFENjV6KqEENWQhLs1aA07v7eMq+dchB4vQLd/gpOr0ZUJIaopeUO1os4mwOxB8PsT4N0cRq+DyJck2IUQhjJNz93m3lAtyIMNMbD6fXB0gds+hg4Py3Z3QghTME2425SkrTBvLKTug9BB0P99qN3A6KqEEOK/JNyvR04mxL4Fm6eDRwMY9j20uM3oqoQQ4m8k3Mvq0GLL0gEXkqHjo9DrVahR2+iqhBCiRBLu15J52rLd3f7fwbcl3P01NLrR6KqEEOKqrB7uSqkg4GXAU2s9pLRjpldYCNu/guVvQH62pafedRw4OhtdmRBCXFOZpnYopWYppVKVUnuLHe+nlDqklIpXSo0H0Fonaq1HFf2+ko6ZWtoh+GqA5S7TBm3hiY3Q/VkJdiGEzSjrvL2vgH5FDyilHIEpQH+gJTBcKdXSqtVVtfwcWDkBPusKqQdg8BR4cD54NTO6MiGEuC5lGpbRWq9RSgUWO9wRiNdaJwIopeYAg4H9Vq2wqhzbYNmc+sxhaHM33DoBavkYXZUQQpRLRe648QdOFHmcBPgrpbyUUtOAMKXUiwAlHSsuLS2NiIiI/35U2Q1Nl89bQv3L/pax9RH/gbtmSLALIWxaRd5QLWntWq21PguMLnbwb8eKc3FxITw8nKioKKKioipQVhlpbZkBs/h5uJQGN42xLBvg4l751xZCiEpWkXBPAhoVeRwAJJf3ZFW6/EBGEix8Fg4vhvpt4d6foGH7qrm2EEJUgYoMy2wBQpRSTZVSLsAwYF55T1YlC4cVFsCf02BKJziyGm55Gx5dKcEuhLA7Zeq5K6V+AHoC3kqpJOB1rfVMpdQYYCngCMzSWu8rbyGV3nM/vRfmj4WT26BZbxj4MdQNrLzrCSGEgZTW2ugaAAgJCdGRkZHWH3PPuwyr34MNk6BGHcvm1G2GyHZ3Qgi7oJTaprWOKH7cNMsPVErPPWGl5Uakc0eg/X1wy1vgVs+61xBCCBMyzeLjVh1zz0qHXx+Hb2639NAfmAe3T5FgF0JUG/bVc9ca9vwMS8ZDdoZlq7vuz4FzTesUKYQQNsI04V5h6Udg4TOQsAL8I2BQDPi1MroqIYQwhO0PyxTkw/qJMPUmOLEZ+n8Ao5ZJsAshqjXT9NzLPSyz5AXYMgOaD4ABH4BngPWLE0IIG2OacC+3To9D0+6WvUxleqMQQgD2MCzjHQwtB0uwCyFEEabpuVfp2jJCCGHnTNNzF0IIYT0S7kIIYYdME+5VsiqkEEJUEzLmLoQQdsg0PfeKsKf/FKQt5mMv7QBpi1lVRltMG+5Fh2eu9XlF/mHKMgx0te8p6WvFj9lCW663HcUf//V50WO22BZrPydXq7Ms32Mvv1+lfc0W22K210ppTBPu3t7e/+/x9TzJFWG2X9iKMNsvbEWYKdwrykzhXhHyWin9uJGvFeBMSQdNs1mHUmoJUDThPYGMMn7uTSkNLIOi5yvP95T0teLHbKEt19uO4o//+rzoMVtsi7Wfk6vVWZbvsZffr9K+ZottMdtr5YzWul/xg6YJdyGEENZjmmEZIYQQ1iPhLoQQdkjCXQgh7JDdh7tSKlQpNU0pNVcp9bjR9VSEUup2pdQXSqnflVK3GF1PeSmlgpRSM5VSc42upTyUUu5Kqa+vPBcjjK6nImz9uSjKjl4f1sksrbVpP4BZQCqwt9jxfsAhIB4YX8ZzOQAz7aQtdY1qi5XbMdfo37HytAu4H4i68vmPRtdujefITM+FFdpi2OvDyu2oUGYZ3uhr/IN0B8KL/oMAjkACEAS4ALuAlkAbYEGxD98rPzMI2ADca+ttufJzHwHhdtAO0wTKdbbrRaD9le/53ujaK9IWMz4XVmiLYa8Pa7XDGpllmrVlSqK1XqOUCix2uCMQr7VOBFBKzQEGa60nAANLOc88YJ5SaiHwfSWWXCprtEUppYB3gcVa6+2VXHKJrPWcmM31tAtIAgKAnZhwaPM627K/isu7LtfTFqXUAQx+fZTmep8Ta2SW6X4xy8AfOFHkcdKVYyVSSvVUSsUopT4HFlV2cdfputoCPAX0AYYopUZXZmHX6XqfEy+l1DQgTCn1YmUXVwGltesX4C6l1GeArSxjWmJbbOi5KKq058Wsr4/SlPacWCWzTN1zL0VJ++mVeieW1noVsKqyiqmg621LDBBTeeWU2/W24yxgCy++Etultb4EPFzVxVRQaW2xleeiqNLaYtbXR2lKa8cqrJBZtthzTwIaFXkcACQbVEtF2Utb7KUdxdlTu6Qt5lOp7bDFcN8ChCilmiqlXIBhwDyDayove2mLvbSjOHtql7TFfCq3HUa/i3yNd5h/AE4BeVj+lxt15fgA4DCWd5pfNrrO6tQWe2mHPbdL2mK+DyPaIQuHCSGEHbLFYRkhhBDXIOEuhBB2SMJdCCHskIS7EELYIQl3IYSwQxLuQghhhyTchRDCDkm4CyGEHZJwF0IIO/R/Is/QRnYO9J8AAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# piro20\n",
    "\n",
    "\n",
    "import astropy.constants as const\n",
    "import numpy as np\n",
    "def planck_lambda(T, Rbb, lamb):\n",
    "    # 一个黑体谱公式\n",
    "    '''\n",
    "    T in the unit of K\n",
    "    Rbb in the unit of Rsun\n",
    "    lamb in the unit of AA\n",
    "    '''\n",
    "    ANGSTROM = 1.0e-8\n",
    "    # convert to cm for planck equation\n",
    "    lamb2 = lamb * ANGSTROM\n",
    "    x = const.h.cgs.value * const.c.cgs.value / (const.k_B.cgs.value * T * lamb2)\n",
    "    x = np.array(x)\n",
    "    Blambda = (2. * const.h.cgs.value * const.c.cgs.value**2 ) /  (lamb2**5. ) / (np.exp(x) - 1. )\n",
    "    # convert back to ANGSTROM   \n",
    "    spec = Blambda*ANGSTROM # in units of erg/cm2/Ang/sr/s\n",
    "    Rbb *= const.R_sun.cgs.value\n",
    "    spec1 = spec * (4. * np.pi * Rbb**2) * np.pi # erg/AA/s\n",
    "    # spec1 *= 1./ (4*np.pi*D**2) to correct for distance\n",
    "    return spec1\n",
    "from scipy.interpolate import interp1d\n",
    "def model_piro20(t_, Menv_=0.1, Renv=200, Eext49 = 10):\n",
    "    \"\"\"\n",
    "    # https://iopscience.iop.org/article/10.3847/1538-4357/abe2b1\n",
    "    公式全部来源于这里，\n",
    "    \"\"\"\n",
    "\n",
    "    \"\"\"\n",
    "    \n",
    "    t_ = np.logspace(-1, 1)\n",
    "    t_ in day\n",
    "    Renv in R_sun\n",
    "    Menv_ in Msun\n",
    "    \"\"\"\n",
    "    # 全部统一为cgs单位\n",
    "    t = t_ * 24*3600 #  in seconds\n",
    "    Ee = 1e+49 * Eext49\n",
    "    Renv= Renv  * const.R_sun.cgs.value\n",
    "    Me  = Menv_ * const.M_sun.cgs.value\n",
    "    # 一些常亮\n",
    "    n = 10\n",
    "    delta = 1.1\n",
    "    K = (n-3) * (3-delta) / (4 * np.pi * (n-delta)) # K = 0.119\n",
    "    kappa = 0.2\n",
    "    c = const.c.cgs.value\n",
    "    \n",
    "    vt = ((n-5) * (5-delta) / ((n-3) * (3-delta)))**0.5 * (2 * Ee / Me)**0.5\n",
    "    td = (3 * kappa * K * Me / ((n-1) * vt * c))**0.5 # in second\n",
    "    #td_day = td / (24*3600)\n",
    "    \n",
    "    prefactor = np.pi*(n-1)/(3*(n-5)) * c * Renv * vt**2 / kappa \n",
    "    L1 = prefactor * (td/t)**(4/(n-2))\n",
    "    L2 = prefactor * np.exp(-0.5 * ((t/td)**2 - 1))\n",
    "    Ls = np.zeros(len(t))\n",
    "    ix1 = t < td\n",
    "    Ls[ix1] = L1[ix1]\n",
    "    Ls[~ix1] = L2[~ix1]\n",
    "    tph = (3 * kappa * K * Me / (2*(n-1)*vt**2))**0.5\n",
    "    R1 = (tph/t)**(2/(n-1)) * vt * t\n",
    "    R2 = ((delta-1)/(n-1)*((t/td)**2-1)+1)**(-1 /(delta+1)) * vt * t\n",
    "    Rs = np.zeros(len(t))\n",
    "    # print((td*u.s).to(u.d))\n",
    "    # print((tph*u.s).to(u.d))\n",
    "    ix1 = t < tph\n",
    "    Rs[ix1] = R1[ix1]\n",
    "    Rs[~ix1] = R2[~ix1]\n",
    "    # if type(wv)==str:\n",
    "    #     from SNeTools.instrument_info import center_wv\n",
    "    #     Luminosity_lambda={}\n",
    "    #     for f in wv:\n",
    "    #         sigmaT4 = Ls / (4 * np.pi * Rs **2)\n",
    "    #         T = (sigmaT4/const.sigma_sb.cgs.value)**(1/4.)\n",
    "    #         Llambda = planck_lambda(T, Rs/const.R_sun.cgs.value, center_wv[f] )\n",
    "    #         ix = t<0\n",
    "    #         Llambda[ix] = 0\n",
    "    #         Luminosity_lambda[f]=Llambda\n",
    "    #         # Rph=Vsc*t\n",
    "    #         # T = ((Luminosity / (4 * np.pi * Rph **2 ) /const.sigma_sb.cgs)**(1/4)).to(u.K)\n",
    "    #         # Luminosity_lambda[i]=planck_lambda(T.value, Rph.cgs/const.R_sun.cgs, center_wv[i])*u.Unit(\"erg AA-1 s-1\")\n",
    "            \n",
    "    #     return Luminosity_lambda\n",
    "\n",
    "    return Ls,Rs,\n",
    "\n",
    "# def photo_radii_p20(t,\n",
    "#                     E  =1e52*u.erg,\n",
    "#                     Mej=1*const.M_sun.cgs,\n",
    "#                     kappa=0.1*u.Unit(\"cm2 g-1\"),\n",
    "#                     # V_t=1e-1*u.Unit(\"km s-1\"),\n",
    "#                     delta=1.1,n=10):\n",
    "#     if not isinstance(E    ,u.Quantity): E  =E*u.erg\n",
    "#     if not isinstance(Mej  ,u.Quantity): Mej=Mej*const.M_sun.cgs\n",
    "#     if not isinstance(kappa,u.Quantity): kappa=kappa*u.Unit(\"cm2 g-1\")\n",
    "#     t=t.cgs\n",
    "#     K   = (n-3)*(3-delta)/(4*np.pi*(n-delta))\n",
    "#     V_t = (((n-5)*(5-delta))/((n-3)*(3-delta)))**(1/2) * ((2*E/Mej)**(1/2)).to(u.Unit('km s-1'))\n",
    "#     t_ph= ((3*kappa*K*Mej)/ (2*(n-1)*V_t**2))**(1/2)\n",
    "    \n",
    "#     r1 = V_t*t * (t_ph/t)**(2/(n-1))\n",
    "#     if delta==1:\n",
    "#         r2= V_t*t * np.e**(1/(n-1))* np.exp(-1/(n-1)*(t/t_ph)**2)\n",
    "#     else:\n",
    "#         r2= V_t*t *(((delta-1)/(n-1))*((t/t_ph)**2-1)+1)**(-1/(delta-1))\n",
    "    \n",
    "#     r = np.zeros(len(t))*u.cm\n",
    "#     t_arg=t.cgs.value<t_ph.cgs.value\n",
    "    \n",
    "#     r[t_arg] =r1[t_arg]\n",
    "#     r[~t_arg]=r2[~t_arg]\n",
    "    \n",
    "#     # rt=t*V_t\n",
    "#     # plt.plot(t.to(u.d),rt.cgs,color=\"b\",label=\"r_br\")         \n",
    "#     # plt.plot(t.to(u.d),r1.cgs,color=\"g\",label=\"r1\")         \n",
    "#     # plt.plot(t.to(u.d),r2.cgs,color=\"r\",label=\"r2\")             \n",
    "#     # plt.plot(t.to(u.d), r.cgs,color=\"k\",label=\"r\",ls=\"--\",)\n",
    "#     # plt.legend()\n",
    "#     return r.cgs    \n",
    "#    \n",
    "def Arnett_NiCo(t, tau_m = 13, Mni = 0.05, t0 = 30, texp=0.0,r=0.99):\n",
    "\n",
    "    import warnings\n",
    "    warnings.filterwarnings(\"ignore\")\n",
    "\n",
    "    # 此部分也可以参考MOSFIT的简介\n",
    "    '''\n",
    "    Calculate the flux of a radioactivity powered SN at photospheric phase\n",
    "    \n",
    "    ts is in the unit of day\n",
    "    Mni_ is in the unit of Msun\n",
    "    \n",
    "    The euqation is from\n",
    "    Valenti 2008 MNRAS 383 1485V, Appendix A\n",
    "    '''\n",
    "    epsilon_Ni56 = 3.9e10*u.Unit(\"erg s-1 g-1\").cgs.scale\n",
    "    epsilon_Co56 = 6.78e9*u.Unit(\"erg s-1 g-1\").cgs.scale\n",
    "    tau_Ni56 = 8.8*u.d.cgs.scale\n",
    "    tau_Co56 = 111.38*u.d.cgs.scale\n",
    "\n",
    "    t     = t     * u.d.cgs.scale if not isinstance(t     ,u.Quantity) else t.cgs.value #  in seconds\n",
    "    texp  = texp  * u.d.cgs.scale if not isinstance(texp  ,u.Quantity) else texp.cgs.value #  in seconds\n",
    "    tau_m = tau_m * u.d.cgs.scale if not isinstance(tau_m ,u.Quantity) else tau_m.cgs.value #  in seconds\n",
    "    t0    = t0    * u.d.cgs.scale if not isinstance(t0    ,u.Quantity) else t0.cgs.value #  in seconds\n",
    "    Mni   = Mni   * const.M_sun.cgs.value if not isinstance(Mni,u.Quantity) else Mni.cgs.value #  in seconds\n",
    "\n",
    "    t=t-texp\n",
    "\n",
    "    L_nico=np.zeros(t.shape)\n",
    "    L_nico=Mni * (\n",
    "            epsilon_Ni56 * np.exp(-t / tau_Ni56) +\n",
    "            epsilon_Co56 * np.exp(-t / tau_Co56))\n",
    "    x=t/tau_m\n",
    "    # 创建积分变量\n",
    "    # 影响积分效率的两个参数\n",
    "    IntNum=1000\n",
    "    MinLogSpace=-3\n",
    "    lsp = np.logspace(-np.log10(max(x)) +MinLogSpace, 0, int(IntNum/2))\n",
    "    xm=np.unique(np.concatenate((lsp, 1 - lsp)))\n",
    "    # print(\"xm\",xm)\n",
    "    int_x = np.clip((x.reshape(len(x),1)) * xm, min(x),max(x))\n",
    "    int_xe = x.reshape(len(t), 1)\n",
    "    # print(min(int_x[:,-1]),max(int_x[:,-1]),)\n",
    "    # print(np.min(int_x))\n",
    "    # print(np.max(int_x))\n",
    "    # print()\n",
    "    int_L=interp1d(int_x[:,-1],L_nico,copy=False,assume_sorted=True)(int_x)\n",
    "    \n",
    "    Linp=int_L*int_x*np.exp(int_x**2-int_xe**2)\n",
    "    l = np.trapz(Linp, int_x)\n",
    "\n",
    "\n",
    "    l=2 * (-np.expm1(-(t0/t)**2))*l\n",
    "\n",
    "    Luminosity=l*(1 + 1.4 * (2 + x/0.59) / (1 +np.exp(x/0.59)) *(r - 1))    \n",
    "\n",
    "    Luminosity[t<0]=0.0\n",
    "    Luminosity[np.isinf(Luminosity)]=0.0\n",
    "    Luminosity[np.isnan(Luminosity)]=0\n",
    "    Luminosity[Luminosity<0        ]=0\n",
    "\n",
    "    return Luminosity# *u.Unit(\"erg s-1\")\n",
    "# t=np.linspace(,200,1000)\n",
    "\n",
    "def blackbody_radiation(T, wavelength):\n",
    "    # 一个黑体谱\n",
    "    h  =const.h\n",
    "    c  =const.c\n",
    "    k_B=const.k_B\n",
    "    # 黑体辐射公式\n",
    "    I = (2 * h * c**2 / wavelength**5) / (np.exp((h * c) / (wavelength * k_B * T)) - 1)\n",
    "\n",
    "    return I.to(u.erg / (u.cm**2 * u.s * u.AA ))/u.sr\n",
    "\n",
    "global center_wv \n",
    "fit_band=[\"o\",\"c\",\"g\",\"r\"]\n",
    "\n",
    "center_wv=np.array([1])*u.AA\n",
    "delta=1.1\n",
    "n=12.5\n",
    "kappa=0.2\n",
    "M_sun=const.M_sun.cgs.value\n",
    "R_sun=const.R_sun.cgs.value\n",
    "def sce_NiCo_multicolor(t,theta,*args):\n",
    "    # 将AL Piro的模型，和Arnett模型加载一起，并计算一个多色光变\n",
    "    Beta=13.8\n",
    "    c=const.c\n",
    "    # 需要拟合的这些参数。\n",
    "    Vphot,\\\n",
    "    texp,\\\n",
    "    tau_m,\\\n",
    "    t_trap,\\\n",
    "    f_Ni,\\\n",
    "    T_floor,\\\n",
    "    Renv,\\\n",
    "    Menv,\\\n",
    "    Eext49=tuple(theta)\n",
    "    t=t-texp\n",
    "    # 抛射物和散射时标的关系。\n",
    "    Mej=(Beta*c*Vphot*u.Unit(\"km s-1\"))*(tau_m*u.d)**2/(2*kappa*u.Unit(\"cm2 g-1\"))\n",
    "    Mej=Mej.cgs.value/M_sun\n",
    "\n",
    "    L_sce,R_sce=model_piro20(t, Menv_=Menv, Renv=Renv, Eext49=Eext49)\n",
    "    L_bol=np.zeros(t.shape)\n",
    "    print(L_sce)\n",
    "    print(Arnett_NiCo(t, tau_m = tau_m, Mni = Mej*f_Ni, t0 = t_trap, texp=0.0,r=0.99))\n",
    "    L_bol+=L_sce\n",
    "    L_bol+=Arnett_NiCo(t, tau_m = tau_m, Mni = Mej*f_Ni, t0 = t_trap, texp=0.0,r=0.99)\n",
    "    \n",
    "    \n",
    "    # 设置光球层\n",
    "    R_ph_Ni=Vphot*1e5*t*86400\n",
    "    R_ph= np.maximum(R_sce, R_sce)\n",
    "\n",
    "    T=(L_bol/(4*np.pi*R_ph**2)/const.sigma_sb.cgs.value)**(1/4)\n",
    "\n",
    "    # K   = (n-3)*(3-delta)/(4*np.pi*(n-delta))\n",
    "    # t_ph= ((3*kappa*K*Mej)/ (2*(n-1)*Vphot**2))**(1/2)\n",
    "    # # print(t_ph.to(u.d))\n",
    "    # # Sucess_LC=True\n",
    "    # # if np.where(L==0)[0].any():Sucess_LC=False\n",
    "    # time_rise_t=t[t.cgs>t_ph.cgs][np.min(np.where(np.diff(T[t.cgs>t_ph.cgs])>0)[0])]\n",
    "    # T[t.cgs>time_rise_t.cgs]=np.nanmin(T[t.cgs<time_rise_t.cgs])\n",
    "    # 假设光球层，用一个T_floor模型\n",
    "    T[T<T_floor]=T_floor\n",
    "    R_ph=np.sqrt(L_bol/(const.sigma_sb*T**4)/(4*np.pi))\n",
    "\n",
    "    # I=blackbody_radiation(T.decompose(),center_wv).value\n",
    "    plt.figure()\n",
    "    plt.plot(t,L_bol)\n",
    "    plt.loglog()\n",
    "    plt.figure()\n",
    "    plt.plot(t,T)\n",
    "    plt.loglog()\n",
    "    plt.figure()\n",
    "    plt.plot(t[1:],np.diff(R_ph)/np.diff(t*86400)/1e5)\n",
    "    plt.loglog()\n",
    "    plt.figure()\n",
    "    plt.plot(t,R_ph)\n",
    "    plt.plot(t,R_ph_Ni)\n",
    "    plt.loglog()\n",
    "    \n",
    "    I=blackbody_radiation(T*u.K,center_wv).value\n",
    "\n",
    "    # I=np.array(I)*u.Unit(\"erg cm-2 s-1 sr-1 AA-1\")\n",
    "    # F=(np.pi*u.sr) * \\\n",
    "    #     I*u.Unit(\"erg cm-2 s-1 sr-1 AA-1\") *\\\n",
    "    #     (rphot*np.sin(theta*deg)/(DL*u.cm))**2\n",
    "    # print(rphot)\n",
    "    # plt.subplot(311)\n",
    "    # plt.semilogy(t,L)\n",
    "    # plt.subplot(312)\n",
    "    # plt.plot(t,rphot)\n",
    "    # plt.subplot(313)\n",
    "    # plt.plot(t,T)\n",
    "    # print(L)\n",
    "    # print(F)\n",
    "    # return F,Sucess_LC\n",
    "    return T\n",
    "\n",
    "t=np.logspace(-3,3,200)\n",
    "# l1=Arnett_NiCo(t,6,0.12,40,-30)\n",
    "# plt.semilogy(t,l1)\n",
    "# plt.loglog()\n",
    "theta=np.array([5000,#Vphot,\\\n",
    "                # 60000,#texp,\\\n",
    "                0.0,#texp,\\\n",
    "                10,#tau_m,\\\n",
    "                30,#t_trap,\\\n",
    "                0.1,#f_Ni,\\\n",
    "                3000,#T_floor,\\\n",
    "                500,#Renv,\\\n",
    "                0.2,#Menv,\\\n",
    "                15])#Eext49\n",
    "\n",
    "_=sce_NiCo_multicolor(t,theta)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# mcmc包\n",
    "# 这个是集成emcee下的一个包。\n",
    "# 如果理解完成上面的东西了之后，可以用这个来拟合数据。\n",
    "# 还没整理好\n",
    "class curve_fit_mcmc():\n",
    "    def __init__(self,f,xdata,ydata,sigma=[],bounds=[],\n",
    "                 p0=[],\n",
    "                 N_pool=4,seed=42,nwalkers=32,\n",
    "                 nsamplers=5000,nburn=1500,save_chain=True,\n",
    "                 prior=None,\n",
    "                 ylimit=np.array([]),args=None\n",
    "                 ):\n",
    "        \"\"\"\n",
    "        f: 拟合的函数\n",
    "        x_data: 数据x\n",
    "        y_data: 数据y\n",
    "        y_err: 数据y的误差\n",
    "        bounds: 参数的边界，默认为空，即不限制参数的取值范围[-inf,inf]\n",
    "        p0: 初始参数\n",
    "        N_pool: 并行计算的进程数\n",
    "        seed: 随机种子\n",
    "        nwalker: 每个walker的参数个数\n",
    "        nsamplers: 采样步数\n",
    "        nburn: 预热步数\n",
    "        save_chain: 是否保存采样结果\n",
    "        ylimit :numpy array. 记录数据是上限(>0)/下限(<0)/实际探测值(0).\n",
    "        prior: 先验方程，规定了函数之间内部的一些关系。\n",
    "\n",
    "        \"\"\"\n",
    "        import inspect\n",
    "        self.func=f\n",
    "        self.x_data=xdata\n",
    "        self.y_data=ydata\n",
    "        self.bounds=bounds\n",
    "        self.y_err =sigma\n",
    "        self.args=args\n",
    "        self.N_pool=N_pool\n",
    "        self.p0 = np.array(p0)\n",
    "        self.seed=seed\n",
    "        self.ndim= len(p0)\n",
    "        self.nwalker=nwalkers\n",
    "        self.nsamplers=nsamplers\n",
    "        self.nburn=nburn\n",
    "        self.save_chain=save_chain\n",
    "        self.ylimit=ylimit if np.any(ylimit) else np.zeros(len(ydata))\n",
    "        self.func_name=self.func.__name__\n",
    "        self.prior=self.prior_func if prior==None else prior\n",
    "        self.check_prior=False if prior==None else True\n",
    "        self.check_limit=True if np.any(np.where(self.ylimit!=0)[0]) else False\n",
    "\n",
    "\n",
    "        self.upperlimit_arg=np.where(self.ylimit>0)[0]\n",
    "        self.lowerlimit_arg=np.where(self.ylimit<0)[0]\n",
    "        self.real_data_arg =np.where(self.ylimit==0)[0]\n",
    "        self.x_data_real=self.x_data\n",
    "        self.y_data_real=self.y_data[self.real_data_arg]\n",
    "        self.y_err_real =self.y_err [self.real_data_arg]\n",
    "        self.check_y_err()\n",
    "        # self.run_mcmc()\n",
    "    def prior_func(self,params): \n",
    "        return 0\n",
    "        \n",
    "    def check_limit_data(self,y_model):\n",
    "        if np.any(self.upperlimit_arg):\n",
    "            if np.any(y_model[self.upperlimit_arg]>self.y_data[self.upperlimit_arg]):\n",
    "                return -np.inf\n",
    "        if np.any(self.lowerlimit_arg):\n",
    "\n",
    "            if np.any(y_model[self.lowerlimit_arg]<self.y_data[self.lowerlimit_arg]):\n",
    "\n",
    "                return -np.inf\n",
    "        return 0\n",
    "    \n",
    "    def check_bound(self):\n",
    "        if not any(list(self.bounds)):\n",
    "            self.bounds=[tuple(self.ndim*[-np.inf]),\n",
    "                         tuple(self.ndim*[ np.inf])]\n",
    "    def check_y_err(self):\n",
    "        self.y_err=np.abs(self.y_err)\n",
    "        if np.any(self.y_err == 0):\n",
    "            self.y_err[self.y_err==0]=1e-5\n",
    "            warnings.warn(\"y_err has 0, and which have set to 1e-5.\")\n",
    "        if np.any(np.isnan(self.y_err)):\n",
    "            self.y_err[np.isnan(self.y_err)]=1\n",
    "            warnings.warn(\"y_err has None, and which have set to 1.\")\n",
    "        if np.any(np.isposinf(self.y_err)):\n",
    "            self.y_err[np.isposinf(self.y_err)]=10\n",
    "            warnings.warn(\"y_err has inf, and which have set to 10.\")\n",
    "\n",
    "    def check_bound(self,params):\n",
    "        for i,p in enumerate(params):\n",
    "            if not self.bounds[0][i]<p<self.bounds[1][i]:\n",
    "                return -np.inf\n",
    "        return 0\n",
    "    def log_prob(self,params,*args):\n",
    "        if np.isneginf(self.check_bound(params)): return -np.inf\n",
    "        if self.check_prior: \n",
    "            \n",
    "            if self.prior (params)==-np.inf:  return -np.inf\n",
    "        y_model=self.func(self.x_data,params,args)\n",
    "        if self.check_limit:\n",
    "            \n",
    "            if np.isneginf(self.check_limit_data(y_model)): \n",
    "                return -np.inf\n",
    "\n",
    "        lh=-np.sum(((self.y_data-y_model)/self.y_err)[self.real_data_arg]**2)\n",
    "        \n",
    "        if np.isnan(lh): \n",
    "            return -np.inf\n",
    "        return lh\n",
    "    def run_mcmc(self):\n",
    "        import emcee\n",
    "        from multiprocessing import Pool\n",
    "        from datetime import datetime\n",
    "        import pickle\n",
    "        np.random.seed(self.seed)\n",
    "        pos =  self.p0+ 1e-3 * np.random.random((self.nwalker,self.ndim))\n",
    "        # pos =  self.p0+ 1e-3 * np.random.randn(self.nwalker,self.ndim)\n",
    "\n",
    "        if self.save_chain:\n",
    "            # 将时间格式化为字符串，精确到分钟\n",
    "            time_str = datetime.now().strftime('%Y-%m-%d %H:%M').replace(\" \",\"\")\n",
    "            self.output_name=filename=f\"{self.func_name}_emcee_output_{time_str}.h5\"\n",
    "            if os.path.exists(filename):\n",
    "                os.remove(filename)\n",
    "            backend = emcee.backends.HDFBackend(filename)\n",
    "        \n",
    "            with Pool(self.N_pool) as pool:\n",
    "                sampler = emcee.EnsembleSampler(self.nwalker, \n",
    "                                                    self.ndim, \n",
    "                                                    self.log_prob, \n",
    "                                                    args=(self.args),\n",
    "                                                    pool=pool,\n",
    "                                                    backend=backend)\n",
    "                _=sampler.run_mcmc(pos,self.nsamplers, progress=True)\n",
    "                print(self.output_name,\"has been saved in\",os.path.abspath(os.path.curdir))\n",
    "\n",
    "\n",
    "        else:\n",
    "            with Pool(self.N_pool) as pool:\n",
    "                sampler = emcee.EnsembleSampler(self.nwalker, \n",
    "                                                    self.ndim, \n",
    "                                                    self.log_prob, \n",
    "                                                    args=(self.args),pool=pool)\n",
    "                _=sampler.run_mcmc(pos,self.nsamplers, progress=True)\n",
    "\n",
    "        self.samplers=sampler\n",
    "\n",
    "    def display_latex_or_text(self,string):\n",
    "        from IPython.display import display, Math\n",
    "        try:\n",
    "            display(Math(string))\n",
    "        except:\n",
    "            print(string)\n",
    "\n",
    "\n",
    "    def display_h5results(h5name=None,ParamterLabel=[],discard=1000,corner_range=0.99):\n",
    "        import emcee\n",
    "        from corner import corner\n",
    "        def display_latex_or_text(string):\n",
    "            from IPython.display import display, Math\n",
    "            try:\n",
    "                display(Math(string))\n",
    "            except:\n",
    "                print(string)\n",
    "\n",
    "        reader = emcee.backends.HDFBackend(h5name)\n",
    "        flat_samples = reader.get_chain(discard=discard, thin=15,flat=True)\n",
    "        ParamterNum  = len(flat_samples[0])\n",
    "        FitResult    = np.array([np.percentile(flat_samples[:,i],[16,50,84]) for i in range(ParamterNum)])\n",
    "        FitResultVal = FitResult[:,1]\n",
    "        FitResultStd = np.diff(FitResult)\n",
    "        if not any(list(ParamterLabel)):\n",
    "            ParamterLabel = [r\"$Param_{%d}$\"%i for i in range(ParamterNum)]\n",
    "\n",
    "        corner(flat_samples,\n",
    "            labels=ParamterLabel,\n",
    "            truths=FitResultVal,\n",
    "            truth_color=\"k\",quantiles=[0.16, 0.5, 0.84], \n",
    "            title_fmt='.2f',show_titles=True,\n",
    "            range=[corner_range]*ParamterNum)\n",
    "        \n",
    "        for i,val in enumerate(FitResultVal):\n",
    "            display_latex_or_text(f\"${ParamterLabel[i]}={val:.2f} + {FitResultStd[i][0]:.2f} - {FitResultStd[i][1]:.2f}$\")\n",
    "        \n",
    "\n",
    "    def display_results(self,ParamterLabel=[],corner_range=0.99):\n",
    "        self.flat_samples = self.samplers.get_chain(discard=self.nburn, thin=15, flat=True)\n",
    "        self.ParamterNum  = len(self.flat_samples[0])\n",
    "        self.FitResult    = np.array([np.percentile(self.flat_samples[:,i],[16,50,84]) for i in range(self.ParamterNum)])\n",
    "        self.FitResultVal = self.FitResult[:,1]\n",
    "        self.FitResultStd = np.diff(self.FitResult)\n",
    "        if not any(list(ParamterLabel)):\n",
    "            import inspect\n",
    "            parameters = inspect.signature(self.func).parameters\n",
    "            ParamterLabel = [param.name for param in parameters.values()][1:]\n",
    "            if len(ParamterLabel)!=self.ParamterNum:\n",
    "                ParamterLabel = [r\"$Param_{%d}$\"%i for i in range(self.ParamterNum)]\n",
    "\n",
    "        self.ParamterLabel=ParamterLabel\n",
    "        from corner import corner\n",
    "        corner(self.flat_samples,\n",
    "            labels=ParamterLabel,\n",
    "            truths=self.FitResultVal,\n",
    "            truth_color=\"k\",quantiles=[0.16, 0.5, 0.84], \n",
    "            title_fmt='.2f',show_titles=True,\n",
    "            range=[corner_range]*self.ParamterNum)\n",
    "        \n",
    "        for i,val in enumerate(self.FitResultVal):\n",
    "            l=ParamterLabel[i].strip(\"$\")\n",
    "            self.display_latex_or_text(f\"${l}={val:.2f} _ {{-{self.FitResultStd[i][0]:.2f}}} ^ {{+{self.FitResultStd[i][1]:.2f}}}$\")\n",
    "\n",
    "    @property\n",
    "    def best_model(self):\n",
    "        return self.func(self.x_data,self.FitResultVal,self.args)\n",
    "    @property\n",
    "    def residual(self):\n",
    "        return self.best_model -self.y_data\n",
    "    @property\n",
    "    def chi2(self):\n",
    "        return np.sum((self.residual/self.y_err)**2)\n",
    "    @property\n",
    "    def reduced_chi2(self):\n",
    "        return np.sum((self.residual/self.y_err)**2)/(len(self.x_data)-self.ndim)\n",
    "    def latex_text(self,ParamterLabel=[]):\n",
    "        if not any(ParamterLabel): ParamterLabel=self.ParamterLabel\n",
    "        for i,val in enumerate(self.FitResultVal):\n",
    "            print(\"%9s = $%.2f_{-%.2f}^{+%.2f}$\"%(ParamterLabel[i],val,self.FitResultStd[i,0],self.FitResultStd[i,1]))\n",
    "\n",
    "def model_piro20_mcmc(t,param,*arg):\n",
    "    Menv, Renv, Eext49 ,n,delta=tuple(param)\n",
    "    return model_piro20(t,Menv, Renv, Eext49 ,n,delta)\n",
    "\n",
    "m=curve_fit_mcmc(f=model_piro20_mcmc,\n",
    "                    xdata=BoloLC_fit.mjd[i]-SN2024ggi_exptime.mjd,\n",
    "                    ydata=BoloLC_fit.loc[i,[\"R\",\"T\"]].values.T.flatten(),\n",
    "                    sigma=BoloLC_fit.loc[i,[\"Re\",\"Te\"]].values.T.flatten(),\n",
    "                    #  p0=[0.1, 11.5, 10,8.7, 1.1],\n",
    "                    #  bounds=([0.01, 10, 1,5,1],[0.9, 14, 1000,15,3])\n",
    "                    p0=[0.1, 11.5, 10,8.7, 1.1],\n",
    "                    bounds=([0.01, 10, 1,5,1],[0.9, 14, 1000,15,3])\n",
    "                 )\n",
    "m.run_mcmc()\n",
    "ParamterLabel=[\n",
    "    r\"$M_{env}$\",\n",
    "    r\"$log_{10}(R_{env})$\",\n",
    "    r\"$E_{inj}$\",\n",
    "    \"n\",\n",
    "    r\"$\\delta$\",\n",
    "   ]\n",
    "\n",
    "m.display_results(ParamterLabel)\n",
    "print(\"Reduced Chi2:\"+\"%.2f\"%m.reduced_chi2)\n",
    "Menv, Renv, Eext49 ,n,delta=tuple(m.FitResultVal)\n",
    "results=model_piro20(t,Menv, Renv, Eext49 ,n,delta)\n",
    "# L_cool=results.reshape(3,-1)[0]\n",
    "R_cool=results.reshape(2,-1)[0]\n",
    "T_cool=results.reshape(2,-1)[1]\n",
    "L_cool=R_cool**2*4*np.pi*T_cool**4*const.sigma_sb.cgs.value\n",
    "print(\"Menv:%.2f Msun\"%p[0])\n",
    "print(\"Renv:%.2f Rsun\"%(10**p[1]/const.R_sun.cgs.value))\n",
    "print(\"Einj:%.2f ×10^49 erg\"%p[2])\n",
    "print(\"n   :%.2f\"%p[3])\n",
    "print(\"delta:%.2f\"%p[4])"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Astro",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
